!##############################################################################
! MODULE simulate_economy
!
! This module simulates the shock paths and economy after 
! solving the life-cycle problem 
!
! Author: Annika Bacher (annika.bacher@bi.no)
!          
! This version: May 2024
!
!##############################################################################

module simulate_economy

        use toolbox_own
        use globals
        
        ! set some variables
        
        implicit none
        
        integer :: tm, zs, zps, tps, m, n, nsize, clock

       integer, allocatable :: seed(:)
        
        real*8, dimension(2) :: coupstate
        
        real*8  :: var
        

contains


subroutine shock_path

        ! NOTE: shock paths are by gender, NOT by marital status

        ! initialize random number draw 
          CALL RANDOM_SEED(size = nsize)
          ALLOCATE(seed(nsize))
          CALL SYSTEM_CLOCK(COUNT=clock)
          
          !! when initializing seed according to current time
          !seed = clock + 37*(/ (i - 1, i = 1, nsize) /)
          
          !! initial seed for repliaction
          seed = (/6944626, 6944663, 6944700, 6944737, 6944774, 6944811, 6944848, 6944885 /)
          call random_seed(put=seed)
  
  
          open(unit=1,file='results/seed.txt', &
          status='replace',action='write')
          write(1,*) seed
          close(unit=i)
		
        ! ****************************************************
        !       AGGREGATE STATE 
        ! ****************************************************
        
         ! state 1: boom, state 2: recession
         
         ! realization of aggregate state differs across individuals
	     ! to avoid aggregate state driving shape of life-cycle profiles

         do tm= 1,TT	
			do s = 1,SS
			
				 call random_number(rr)
                 if(rr<=p_rec) then
                 chain_agg(s,tm) = 2     
				 else
				 chain_agg(s,tm) = 1
				 end if
	
			 end do
		 end do 
           
        ! ****************************************************
        !       MARITAL STATUS AND PRODUCTIVITY TRANSITION
        ! **************************************************** 
            
        ! initial productivity: 
		! 50% in lowest persistent, 50% in medium persistent
     
        do s = 1,SS
         
            !  labor productivity shocks female
            call random_number(rr)

             if(rr<0.5d0) then
             z0f(s)=2 
             else 
             z0f(s)=5
             end if 
            
            !  labor productivity shocks male
            call random_number(rr)            

			 if(rr<0.5d0) then
             z0m(s)=2
             else 
             z0m(s)=5 
             end if
            
        end do
            
        do s = 1,2*SS
            ! labor productivity shocks couples
            call random_number(rr)

		     if(rr<0.5d0) then
             z0c(s)=8
             else 
             z0c(s)=36
             end if

        end do 
          
      !!!!!!!!!!!!!!!!!!! chain of shocks !!!!!!!!!!!!!!!!!!!!
         
            ! labor productivity shock 
            
            ! first period -- single women
            do s = 1, SS
            chain_prodf(s,1) = z0f(s)
            end do
            
            ! first period -- single men
            do s = 1, SS
            chain_prodm(s,1) = z0m(s)
            end do
            
            ! first period -- couples
            do s = 1, 2*SS
            chain_prodc(s,1) = z0c(s)
            end do
            
            ! for couples: "translating" into individual state variable  
            do s = 1, SS
            
            do m=1,nz
            do j= 1,nzz
            do i= 1,nz
            
            if(chain_prodc(s,1) == (m-1)*nzz*nz+(j-1)*nz+i) then
            chain_prodmf(s,1) = i +(m-1)*nz
            chain_prodmm(s,1) = j
            end if
            
            end do 
            end do
            end do

             end do
                
            !! marital transitions !!
            
            ! seperate low/ high ability 
            ! women
            AB2_f = dist_thetaf(1)*SS 
            AB_f = nint(AB2_f)
            
            ! men 
            AB2_m = dist_thetam(1)*SS 
            AB_m = nint(AB2_m)


            ! initial period -- marital status (if and to whom married)
            ! uniform distribution across partners
                
                 do s = 1, AB_f    ! low ability women
                 call random_number(rr)
                  
                 if(rr<=0.3592d0) then
                 chain_martf(s,1) = 0                               ! single
                 else if(rr>0.3592d0 .AND. rr<=0.5322d0) then
                 call random_number(aa)
                 do i = 1,9
                 if(aa>((i-1)*(1d0/9d0))) chain_martf(s,1) = i+9    ! married to high ability
                 end do
                 else if (rr>0.5322d0 .AND. rr<=1d0) then
                 call random_number(aa)
                 do i = 1,9
                 if(aa>((i-1)*(1d0/9d0))) chain_martf(s,1) = i      ! married to low ability
                 end do
                 end if 
                
                 end do 
                 
                 do s = AB_f+1, SS ! high ability women
                 call random_number(rr)
                  
                 if(rr<=0.3574d0) then
                 chain_martf(s,1) = 0                              ! single
                 else if(rr>0.3574d0 .AND. rr<=0.5245d0) then
                 call random_number(aa)
                 do i = 1,9
                 if(aa>((i-1)*(1d0/9d0))) chain_martf(s,1) = i     ! married to low ability
                 end do
                 else if (rr>0.5245d0 .AND. rr<=1d0) then
                 call random_number(aa)
                 do i = 1,9
                 if(aa>((i-1)*(1d0/9d0))) chain_martf(s,1) = i+9   ! married to high ability
                 end do
                 end if        
                 
                 end do
                 

                 do s = 1, AB_m ! low ability men  
                 call random_number(rr)
                
                 if(rr<=0.3382) then
                 chain_martm(s,1) = 0                              ! single
                 else if(rr>0.3382 .AND. rr<=0.7419d0) then
                 call random_number(aa)
                 do i = 1,9
                 if(aa>((i-1)*(1d0/9d0))) chain_martm(s,1) = i     ! married to low ability
                 end do
                 else if (rr>0.7419d0 .AND. rr<=1d0) then
                 call random_number(aa)
                 do i = 1,9
                 if(aa>((i-1)*(1d0/9d0))) chain_martm(s,1) = i+9   ! married to high ability
                 end do
                 end if 
                 
                 end do
                 
                 do s = AB_m+1, SS ! high ability men
                 
                 call random_number(rr)
                 if(rr<=0.3218d0) then
                 chain_martm(s,1) = 0                              ! single
                 else if(rr>0.3218d0 .AND. rr<=0.8914d0) then
                 call random_number(aa)
                  do i = 1,9
                 if(aa>((i-1)*(1d0/9d0))) chain_martm(s,1) = i+9   ! married to high ability
                 end do
                 else if (rr>0.8914d0 .AND. rr<=1d0) then
                 call random_number(aa)
                  do i = 1,9
                 if(aa>((i-1)*(1d0/9d0))) chain_martm(s,1) = i     ! married to low ability
                 end do
                 end if
                 
                 end do
      
        write(*,*) 'first period of income and marital shock path done' 
            
                
        !!!!!!!!!!!!! all consecutive periods !!!!!!!!!!!!!
            
            
        ! ***********************
        !       ASSET RETURN
        ! ***********************
            
            do s = 1,SS
				do tm = 1,TT
            
            if (chain_agg(s,tm) == 1) then    ! boom 
				call random_number(rr)
				ZC_st = cumsum(pir_boom(:),nrr) 
				do j=nrr,1,-1
                if(rr<ZC_st(j)) chain_retf(s,tm)=j 
				end do
				
				call random_number(rr)
				ZC_st = cumsum(pir_boom(:),nrr) 
				do j=nrr,1,-1
                if(rr<ZC_st(j)) chain_retm(s,tm)=j 
				end do
				
			else                             ! recession
				call random_number(rr)
				ZC_st = cumsum(pir_rec(:),nrr) 
				do j=nrr,1,-1
                if(rr<ZC_st(j)) chain_retf(s,tm)=j 
				end do
				
				call random_number(rr)
				ZC_st = cumsum(pir_rec(:),nrr) 
				do j=nrr,1,-1
                if(rr<ZC_st(j)) chain_retm(s,tm)=j 
				end do
			end if
			
				end do
			end do
            


		! **********************************************
        !       MARITAL AND PRODUCITIVITY CHAINS 
        ! **********************************************

        ! ********************* WOMEN ******************
            
            !!!!!!!!!  low ability type  !!!!!!!!
            
             do s = 1, AB_f
                do tm = 2,TT
             
            if(chain_martf(s,tm-1) == 0) then            ! LAST PERIOD: single      

            call random_number(rr)
            if(rr<mu_f(tm-1,1,chain_prodf(s,tm-1))) then ! getting married, determine to which type acc. to trans_part
            ZC_mf = cumsum(trans_part_f(chain_prodf(s,tm-1),:),naa*nzz) 
            call random_number(rr) 
            do i=naa*nzz,1,-1
            if(rr<ZC_mf(i)) chain_martf(s,tm) = i 
            end do 
            
            ! productivity 
            if (chain_agg(s,tm) == 1) then    ! boom 
            call random_number(rr)
            ZC_zf = cumsum(pifs_boom(chain_prodf(s,tm-1),:),nzz) 
            do i=nzz,1,-1
                if(rr<ZC_zf(i)) chain_prodmf(s,tm)=i
            end do 
            else 							  ! recession
            call random_number(rr)
            ZC_zf = cumsum(pifs_rec(chain_prodf(s,tm-1),:),nzz) 
            do i=nzz,1,-1
                if(rr<ZC_zf(i)) chain_prodmf(s,tm)=i
            end do 
            end if 
            
            ! translating into chain_prodc 
            do i = 1,nzz
                do j=1,nzz
            if 	   (chain_prodmf(s,tm)== i .AND.  chain_martf(s,tm) == j .AND. i < 4 .OR. chain_prodmf(s,tm)== i .AND.  chain_martf(s,tm) == j+9 .AND. i < 4) then
                chain_prodc(s,tm) = (j-1)*nz+i
            elseif (chain_prodmf(s,tm)== i .AND.  chain_martf(s,tm) == j .AND. i >= 4 .AND. i <= 6 .OR. chain_prodmf(s,tm)== i .AND.  chain_martf(s,tm) == j+9 .AND. i >= 4 .AND. i <= 6) then
                chain_prodc(s,tm) = (nzz*nz)+(j-1)*nz+i-nz
            elseif (chain_prodmf(s,tm)== i .AND.  chain_martf(s,tm) == j .AND. i > 6 .OR. chain_prodmf(s,tm)== i .AND. chain_martf(s,tm) == j+9 .AND. i > 6) then
                chain_prodc(s,tm) = 2*(nzz*nz)+(j-1)*nz+i-2*nz
            end if
                end do
            end do
            
            else                              ! staying single
            chain_martf(s,tm) = 0
           
         
            ! productivity 
            if (chain_agg(s,tm) == 1) then    ! boom
            call random_number(rr)
            ZC_zf = cumsum(pifs_boom(chain_prodf(s,tm-1),:),nzz) 
            do i=nzz,1,-1
                if(rr<ZC_zf(i)) chain_prodf(s,tm)=i
            end do 
			else 							  ! recession 
              call random_number(rr)
            ZC_zf = cumsum(pifs_rec(chain_prodf(s,tm-1),:),nzz) 
            do i=nzz,1,-1
                if(rr<ZC_zf(i)) chain_prodf(s,tm)=i
            end do 
            end if
            
            
            end if ! end if statement whether getting married or staying single
         
         
            else                                          ! LAST PERIOD: married                          
            
            do i = 1,9                                    ! LAST PERIOD: married to low ability type
            if (chain_martf(s,tm-1) == i) then
            call random_number(rr)
            if(rr<lambda(tm-1,1,1,chain_prodmf(s,tm-1),i)) then    ! getting divorced 
            chain_martf(s,tm) = 0
            
            ! productivity 
            if (chain_agg(s,tm) == 1) then    ! boom
            call random_number(rr)
            ZC_zc = cumsum(pic_boom(chain_prodc(s,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s,tm)=j 
            end do
            else 							  ! recession 
            call random_number(rr)
            ZC_zc = cumsum(pic_rec(chain_prodc(s,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s,tm)=j 
            end do
            end if
            
            
            ! "translating" into female state variable
            do m=1,nz
            do j= 1,nzz
            do k= 1,nz
            
            if(chain_prodc(s,tm) == (m-1)*nzz*nz+(j-1)*nz+k) then
            chain_prodf(s,tm) = k +(m-1)*nz
            end if
            
            end do 
            end do
            end do

            
            else                    ! staying married
                                    ! both spouses productivities evolve acc. to pic
            
            ! productivity 
            if (chain_agg(s,tm) == 1) then    ! boom
            call random_number(rr)
            ZC_zc = cumsum(pic_boom(chain_prodc(s,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s,tm)=j 
            end do
            else							  ! recession
             call random_number(rr)
            ZC_zc = cumsum(pic_rec(chain_prodc(s,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s,tm)=j 
            end do
            end if
            
            
            ! "translating" into female prod and marital variable (i.e married to which type)       
            do m=1,nz
            do j= 1,nzz
            do k= 1,nz
            
            if(chain_prodc(s,tm) == (m-1)*nzz*nz+(j-1)*nz+k) then
            chain_prodmf(s,tm) = k +(m-1)*nz
            chain_martf(s,tm) = j
            end if
            
            end do 
            end do
            end do
           
            
            end if               ! end if statement whether getting divorced or not  
            end if               ! end if statement if married to low ability type (last period)
            
            end do               ! end do statement for looping over i=1,9
            
            
            
            do i = 10,18                                ! LAST PERIOD: married to high ability type
            if (chain_martf(s,tm-1) == i) then
            call random_number(rr)
            if(rr<lambda(tm-1,1,2,chain_prodmf(s,tm-1),i-9)) then   ! getting divorced 
            chain_martf(s,tm) = 0
            
            ! productivity 
            if (chain_agg(s,tm) == 1) then    ! boom
            call random_number(rr)
            ZC_zc = cumsum(pic_boom(chain_prodc(s,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s,tm)=j 
            end do
            else 							  ! recession
            call random_number(rr)
            ZC_zc = cumsum(pic_rec(chain_prodc(s,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s,tm)=j 
            end do
            end if
            
            
            ! "translating" into female state variable
            do m=1,nz
            do j= 1,nzz
            do k= 1,nz
            
            if(chain_prodc(s,tm) == (m-1)*nzz*nz+(j-1)*nz+k) then
            chain_prodf(s,tm) = k +(m-1)*nz
            end if
            
            end do 
            end do
            end do
            
            else                                                   ! staying married
           
            ! productivity 
            if (chain_agg(s,tm) == 1) then    ! boom
            call random_number(rr)
            ZC_zc = cumsum(pic_boom(chain_prodc(s,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s,tm)=j 
            end do
            else 							  ! recession
            call random_number(rr)
            ZC_zc = cumsum(pic_rec(chain_prodc(s,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s,tm)=j 
            end do
            end if
            
            ! "translating" into female state variable and which type married to
            do m=1,nz
            do j= 1,nzz
            do k= 1,nz
            
            if(chain_prodc(s,tm) == (m-1)*nzz*nz+(j-1)*nz+k) then
            chain_prodmf(s,tm) = k +(m-1)*nz
            chain_martf(s,tm)  = j+9
            end if
            
            end do 
            end do
            end do

    
            end if
            end if
            end do
            
            
            end if  ! end if statement whether single or married last period
            
            end do  ! end loop over individuals (SS)
                end do  ! end loop time (TT)
                
                
            
            !!!!!!!!!!!!!!!!!!!!  women, high ability type  !!!!!!!!!!!!!!!!!!!!
            
             do s = AB_f+1, SS
                do tm = 2,TT
    
            if(chain_martf(s,tm-1) == 0) then            ! LAST PERIOD: single
            call random_number(rr)
            if(rr<mu_f(tm-1,2,chain_prodf(s,tm-1))) then  ! getting married, determine to which type acc. to trans_part
            ZC_mf = cumsum(trans_part_f(nzz+chain_prodf(s,tm-1),:),naa*nzz)
            call random_number(rr) 
            do i=naa*nzz,1,-1
            if(rr<ZC_mf(i)) chain_martf(s,tm) = i
            end do 
           
            ! productivity 
            if (chain_agg(s,tm) == 1) then    ! boom
            call random_number(rr)
            ZC_zf = cumsum(pifs_boom(chain_prodf(s,tm-1),:),nzz) 
            do i=nzz,1,-1
                if(rr<ZC_zf(i)) chain_prodmf(s,tm)=i
            end do 
            else 							  ! recession
             call random_number(rr)
            ZC_zf = cumsum(pifs_rec(chain_prodf(s,tm-1),:),nzz) 
            do i=nzz,1,-1
                if(rr<ZC_zf(i)) chain_prodmf(s,tm)=i
            end do 
            end if 
            
           ! translating into chain_prodc 
            do i = 1,nzz
                do j=1,nzz
            if 	   (chain_prodmf(s,tm)== i .AND. chain_martf(s,tm) == j .AND. i < 4 .OR. chain_prodmf(s,tm)== i .AND. chain_martf(s,tm) == j+9 .AND. i < 4 ) then
                chain_prodc(s,tm) = (j-1)*nz+i
            elseif (chain_prodmf(s,tm)== i .AND. chain_martf(s,tm) == j .AND. i >= 4 .AND. i <= 6 .OR. chain_prodmf(s,tm)== i .AND.  chain_martf(s,tm) == j+9 .AND. i >= 4 .AND. i <= 6) then
                chain_prodc(s,tm) = (nzz*nz)+(j-1)*nz+i-nz
            elseif (chain_prodmf(s,tm)== i .AND. chain_martf(s,tm) == j .AND. i > 6 .OR. chain_prodmf(s,tm)== i .AND. chain_martf(s,tm) == j+9 .AND. i >  6) then
                chain_prodc(s,tm) = 2*(nzz*nz)+(j-1)*nz+i-2*nz
            end if
                end do
            end do

            else
            chain_martf(s,tm) = 0                         ! staying single
            
            ! productivity 
            if (chain_agg(s,tm) == 1) then    ! boom
            call random_number(rr)
            ZC_zf = cumsum(pifs_boom(chain_prodf(s,tm-1),:),nzz) 
            do i=nzz,1,-1
                if(rr<ZC_zf(i)) chain_prodf(s,tm)=i
            end do 
            else 							  ! recession
            call random_number(rr)
            ZC_zf = cumsum(pifs_rec(chain_prodf(s,tm-1),:),nzz) 
            do i=nzz,1,-1
                if(rr<ZC_zf(i)) chain_prodf(s,tm)=i
            end do  
            end if
             
            end if  ! end if statement whether getting married or staying single
            
            
            else                                          ! LAST PERIOD: married  
            
            do i = 1,9                                    ! LAST PERIOD: married to low ability type
            if (chain_martf(s,tm-1) == i) then
            call random_number(rr)
            if(rr<lambda(tm-1,2,1,chain_prodmf(s,tm-1),i)) then  ! getting divorced 
            chain_martf(s,tm) = 0 
            
            ! productivity 
            if (chain_agg(s,tm) == 1) then    ! boom
            call random_number(rr)
            ZC_zc = cumsum(pic_boom(chain_prodc(s,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s,tm)=j 
            end do
            else 							  ! recession
             call random_number(rr)
            ZC_zc = cumsum(pic_rec(chain_prodc(s,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s,tm)=j 
            end do
            end if
            
            ! "translating" into female state variable
            do m=1,nz
            do j= 1,nzz
            do k= 1,nz
            
            if(chain_prodc(s,tm) == (m-1)*nzz*nz+(j-1)*nz+k) then
            chain_prodf(s,tm) = k +(m-1)*nz
            end if
            
            end do 
            end do
            end do
            

            else                                                ! staying married
            
            ! productivity 
            if (chain_agg(s,tm) == 1) then    ! boom
            call random_number(rr)
            ZC_zc = cumsum(pic_boom(chain_prodc(s,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s,tm)=j 
            end do
            else 							  ! recession
            call random_number(rr)
            ZC_zc = cumsum(pic_rec(chain_prodc(s,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s,tm)=j 
            end do
			end if 
            
                  
            ! "translating" into female prod and marital variable (i.e married to which type)       
            do m=1,nz
            do j= 1,nzz
            do k= 1,nz
            
            if(chain_prodc(s,tm) == (m-1)*nzz*nz+(j-1)*nz+k) then
            chain_prodmf(s,tm) = k +(m-1)*nz
            chain_martf(s,tm) = j
            end if
            
            end do 
            end do
            end do
            

            end if
            end if
            
            end do
            
            do i = 10,18                                   ! LAST PERIOD: married to high ability type
            if (chain_martf(s,tm-1) == i) then
            call random_number(rr)
            if(rr<lambda(tm-1,2,2,chain_prodmf(s,tm-1),i-9)) then  ! getting divorced
            chain_martf(s,tm) = 0
            
            ! productivity 
            if (chain_agg(s,tm) == 1) then    ! boom
            call random_number(rr)
            ZC_zc = cumsum(pic_boom(chain_prodc(s,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s,tm)=j
            end do
            else 							  ! recession
             call random_number(rr)
            ZC_zc = cumsum(pic_rec(chain_prodc(s,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s,tm)=j
            end do
            end if
            
            ! "translating" into female state variable
            do m=1,nz
            do j= 1,nzz
            do k= 1,nz
            
            if(chain_prodc(s,tm) == (m-1)*nzz*nz+(j-1)*nz+k) then
            chain_prodf(s,tm) = k +(m-1)*nz
            end if
            
            end do 
            end do
            end do 

            else     ! staying married
            
            ! productivity 
            if (chain_agg(s,tm) == 1) then    ! boom
            call random_number(rr)
            ZC_zc = cumsum(pic_boom(chain_prodc(s,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s,tm)=j 
            end do
            else 							  ! recession
             call random_number(rr)
            ZC_zc = cumsum(pic_rec(chain_prodc(s,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s,tm)=j 
            end do
            end if 
             
             
            ! "translating" into female state variable and which type married to
            do m=1,nz
            do j= 1,nzz
            do k= 1,nz
            
            if(chain_prodc(s,tm) == (m-1)*nzz*nz+(j-1)*nz+k) then
            chain_prodmf(s,tm) = k +(m-1)*nz
            chain_martf(s,tm)  = j+9
            end if
            
            end do 
            end do
            end do
    
            
            end if
            end if
            end do
            
            end if          ! end if statement whether single or married last period
            
            end do          ! end loop over individuals (SS)
                end do      ! end loop time (TT)
                
                
            !		  ***********************************************    
            !         ********************* MEN *********************
            !		  ***********************************************
            
            !!!!!!!!!!!!! low ability type  !!!!!!!!!!!!!!!!
        
           do s = 1, AB_m
                do tm = 2,TT
                
            if(chain_martm(s,tm-1) == 0) then               ! LAST PERIOD: single
            call random_number(rr)
            if(rr<mu_m(tm-1,1,chain_prodm(s,tm-1))) then    ! getting married, determine to which type acc. to trans_part
            ZC_mm = cumsum(trans_part_m(chain_prodm(s,tm-1),:),naa*nzz)
            call random_number(rr) 
            do i=naa*nzz,1,-1
            if(rr<ZC_mm(i)) chain_martm(s,tm) = i 
            end do 
            
            ! productivity 
            if (chain_agg(s,tm) == 1) then    ! boom
            call random_number(rr)
            ZC_zm = cumsum(pims_boom(chain_prodm(s,tm-1),:),nzz) 
            do i=nzz,1,-1
                if(rr<ZC_zm(i)) chain_prodmm(s,tm)=i
            end do 
            else 							  ! recession
            call random_number(rr)
            ZC_zm = cumsum(pims_rec(chain_prodm(s,tm-1),:),nzz) 
            do i=nzz,1,-1
                if(rr<ZC_zm(i)) chain_prodmm(s,tm)=i
            end do 
            end if  
            
            ! translating into chain_prodc
            do i = 1,nzz
                do j=1,nzz
            if 	   (chain_martm(s,tm)== i .AND.  chain_prodmm(s,tm) == j .AND. i < 4 .OR. chain_martm(s,tm)== i+9 .AND. chain_prodmm(s,tm) == j .AND. i < 4) then
                chain_prodc(s+SS,tm) = (j-1)*nz+i
            elseif (chain_martm(s,tm)== i .AND.  chain_prodmm(s,tm) == j .AND. i >= 4 .AND. i <= 6 .OR. chain_martm(s,tm)== i+9 .AND.  chain_prodmm(s,tm) == j .AND. i >= 4 .AND. i <= 6 ) then
                chain_prodc(s+SS,tm) = (nzz*nz)+(j-1)*nz+i-nz
            elseif (chain_martm(s,tm)== i .AND.  chain_prodmm(s,tm) == j .AND. i > 6 .OR. chain_martm(s,tm)== i+9 .AND.  chain_prodmm(s,tm) == j .AND. i >  6 ) then
                chain_prodc(s+SS,tm) = 2*(nzz*nz)+(j-1)*nz+i-2*nz
            end if
                end do
            end do
                  

            else                                            ! staying single
            chain_martm(s,tm) = 0
                        
            !productivity 
            if (chain_agg(s,tm) == 1) then    ! boom
            call random_number(rr)
            ZC_zm = cumsum(pims_boom(chain_prodm(s,tm-1),:),nzz) 
            do i=nzz,1,-1
                if(rr<ZC_zm(i)) chain_prodm(s,tm)=i
            end do 
            else 							  ! recession
			call random_number(rr)
            ZC_zm = cumsum(pims_rec(chain_prodm(s,tm-1),:),nzz) 
            do i=nzz,1,-1
                if(rr<ZC_zm(i)) chain_prodm(s,tm)=i
            end do 
            end if 
            
            end if                         ! end if statement whether getting married or staying single 
            
            else                                            ! LAST PERIOD: married
            
            
            do i = 1,9                                      ! LAST PERIOD: married to low ability type
            if (chain_martm(s,tm-1) == i) then
            call random_number(rr)
            if(rr<lambda(tm-1,1,1,i,chain_prodmm(s,tm-1))) then  ! getting divorced 
            chain_martm(s,tm) = 0
            
            ! productivity 
            if (chain_agg(s,tm) == 1) then    ! boom
            call random_number(rr)
            ZC_zc = cumsum(pic_boom(chain_prodc(s+SS,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s+SS,tm)=j 
            end do
            else 							  ! recession
            call random_number(rr)
            ZC_zc = cumsum(pic_rec(chain_prodc(s+SS,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s+SS,tm)=j 
            end do
            end if
            
            
            ! "translating" into male state variable
            do m=1,nz
            do j= 1,nzz
            do k= 1,nz
            
            if(chain_prodc(s+SS,tm) == (m-1)*nzz*nz+(j-1)*nz+k) then
            chain_prodm(s,tm) = j
            end if
            
            end do 
            end do
            end do
            
            
            else    ! staying married
            
            
            ! productivity 
            if (chain_agg(s,tm) == 1) then    ! boom
            call random_number(rr)
            ZC_zc = cumsum(pic_boom(chain_prodc(s+SS,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s+SS,tm)=j 
            end do
            else 							  ! recession
             call random_number(rr)
            ZC_zc = cumsum(pic_rec(chain_prodc(s+SS,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s+SS,tm)=j 
            end do
            end if 
            
            ! "translating" into male state variable and which type married to
            do m=1,nz
            do j= 1,nzz
            do k= 1,nz
            
            if(chain_prodc(s+SS,tm) == (m-1)*nzz*nz+(j-1)*nz+k) then
            chain_prodmm(s,tm) = j
            chain_martm(s,tm) = k +(m-1)*nz
            end if
            
            end do 
            end do
            end do  
            
            end if                      ! end if statement whether getting divorced or not  
            end if                      ! end if statement if married to low ability type (last period)
            
            end do                      ! end do statement for looping over i=1,9
            
            
            do i = 10,18                                    ! LAST PERIOD: married to high ability type
            if (chain_martm(s,tm-1) == i) then
            call random_number(rr)
            if(rr<lambda(tm-1,2,1,i-9,chain_prodmm(s,tm-1))) then  ! getting divorced 
            chain_martm(s,tm) = 0
            
            ! productivity 
            if (chain_agg(s,tm) == 1) then    ! boom
            ZC_zc = cumsum(pic_boom(chain_prodc(s+SS,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s+SS,tm)=j 
            end do
            else 							  ! recession
             call random_number(rr)
            ZC_zc = cumsum(pic_rec(chain_prodc(s+SS,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s+SS,tm)=j 
            end do
            end if
            
            ! "translating" into male state variable
            do m=1,nz
            do j= 1,nzz
            do k= 1,nz
            
            if(chain_prodc(s+SS,tm) == (m-1)*nzz*nz+(j-1)*nz+k) then
            chain_prodm(s,tm) = j
            end if
            
            end do 
            end do
            end do
            
            else   ! staying married
           
            ! productivity
            if (chain_agg(s,tm) == 1) then    ! boom 
            call random_number(rr)
            ZC_zc = cumsum(pic_boom(chain_prodc(s+SS,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s+SS,tm)=j 
            end do
            else 							  ! recession
            call random_number(rr)
            ZC_zc = cumsum(pic_rec(chain_prodc(s+SS,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s+SS,tm)=j 
            end do
            end if
            
            
            ! "translating" into male state variable
            do m=1,nz
            do j= 1,nzz
            do k= 1,nz
            
            if(chain_prodc(s+SS,tm) == (m-1)*nzz*nz+(j-1)*nz+k) then
            chain_prodmm(s,tm) = j
            chain_martm(s,tm)  = (k +(m-1)*nz) + 9
            end if
            
            end do 
            end do
            end do 
             

            end if
            end if
            end do
            
            end if    ! end if statement whether single or married last period
            
            end do
                end do
  

            !!!!!!!!!!!!!!!!!!!!  men, high ability type  !!!!!!!!!!!!!!!!!!!!

            do s = AB_m+1, SS
             do tm = 2,TT
                
            if(chain_martm(s,tm-1) == 0) then            ! LAST PERIOD: single
            call random_number(rr)
            if(rr<mu_m(tm-1,2,chain_prodm(s,tm-1))) then ! getting married, determine to which type acc. to trans_part
            ZC_mm = cumsum(trans_part_m(nzz+chain_prodm(s,tm-1),:),naa*nzz) 
            call random_number(rr) 
            do i=naa*nzz,1,-1
            if(rr<ZC_mm(i)) chain_martm(s,tm) = i 
            end do 
            
            !productivity 
            if (chain_agg(s,tm) == 1) then    ! boom 
            call random_number(rr)
            ZC_zm = cumsum(pims_boom(chain_prodm(s,tm-1),:),nzz) 
            do i=nzz,1,-1
                if(rr<ZC_zm(i)) chain_prodmm(s,tm)=i
            end do 
            else 							  ! recession
            call random_number(rr)
            ZC_zm = cumsum(pims_rec(chain_prodm(s,tm-1),:),nzz) 
            do i=nzz,1,-1
                if(rr<ZC_zm(i)) chain_prodmm(s,tm)=i
            end do 
            end if 
                  
            ! translating into chain_prodc 
            do i = 1,nzz
                do j=1,nzz
            if 	   (chain_martm(s,tm)== i .AND. chain_prodmm(s,tm) == j .AND. i < 4 .OR. chain_martm(s,tm)== i+9 .AND.  chain_prodmm(s,tm) == j .AND. i < 4) then
                chain_prodc(s+SS,tm) = (j-1)*nz+i
            elseif (chain_martm(s,tm)== i .AND. chain_prodmm(s,tm) == j .AND. i >= 4 .AND. i <= 6 .OR.  chain_martm(s,tm)== i+9 .AND.  chain_prodmm(s,tm) == j .AND. i >= 4 .AND. i <= 6 ) then
                chain_prodc(s+SS,tm) = (nzz*nz)+(j-1)*nz+i-nz
            elseif (chain_martm(s,tm)== i .AND. chain_prodmm(s,tm) == j .AND. i >  6 .OR. chain_martm(s,tm)== i+9 .AND.  chain_prodmm(s,tm) == j .AND. i >  6) then
                chain_prodc(s+SS,tm) = (nzz*nz)+(j-1)*nz+i-2*nz
            end if
                end do
            end do 

            else                                          ! staying single
            chain_martm(s,tm) = 0
            
            !productivity
            if (chain_agg(s,tm) == 1) then    ! boom 
            call random_number(rr)
            ZC_zm = cumsum(pims_boom(chain_prodm(s,tm-1),:),nzz) 
            do i=nzz,1,-1
            if(rr<ZC_zm(i)) chain_prodm(s,tm)=i
            end do 
            else 							  ! recession
            call random_number(rr)
            ZC_zm = cumsum(pims_rec(chain_prodm(s,tm-1),:),nzz) 
            do i=nzz,1,-1
            if(rr<ZC_zm(i)) chain_prodm(s,tm)=i
            end do 
            end if 

            end if						! end if statement whether getting married or staying single

            else                                        ! LAST PERIOD: married
            
            
            do i = 1,9                                  ! LAST PERIOD: married to low ability type
            if (chain_martm(s,tm-1) == i) then
             call random_number(rr)
            if(rr<lambda(tm-1,1,2,i,chain_prodmm(s,tm-1))) then  ! getting divorced 
            chain_martm(s,tm) = 0
            
            ! productivity 
            if (chain_agg(s,tm) == 1) then    ! boom
            call random_number(rr)
            ZC_zc = cumsum(pic_boom(chain_prodc(s+SS,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s+SS,tm)=j 
            end do
            else 							  ! recession
            call random_number(rr)
            ZC_zc = cumsum(pic_rec(chain_prodc(s+SS,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s+SS,tm)=j 
            end do
            end if
            
            ! "translating" into male state variable
            do m=1,nz
            do j= 1,nzz
            do k= 1,nz
            
            if(chain_prodc(s+SS,tm) == (m-1)*nzz*nz+(j-1)*nz+k) then
            chain_prodm(s,tm) = j
            end if
            
            end do 
            end do
            end do   
             
                         
            else                                                ! staying married
            
            ! productivity
            if (chain_agg(s,tm) == 1) then    ! boom 
            call random_number(rr)
            ZC_zc = cumsum(pic_boom(chain_prodc(s+SS,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s+SS,tm)=j 
            end do
            else 							  ! recession
            call random_number(rr)
            ZC_zc = cumsum(pic_rec(chain_prodc(s+SS,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s+SS,tm)=j 
            end do
            end if
            
            ! "translating" into male state variable
            do m=1,nz
            do j= 1,nzz
            do k= 1,nz
            
            if(chain_prodc(s+SS,tm) == (m-1)*nzz*nz+(j-1)*nz+k) then
            chain_prodmm(s,tm) = j
            chain_martm(s,tm) = k +(m-1)*nz
            end if
            
            end do 
            end do
            end do    
            
        
            end if          ! end if statement whether getting divorced or not 
            end if          ! end if statement if married to low ability type (last period)
            
            end do          ! end do statement for looping over i=1,9
            
            
            do i = 10,18                                 ! LAST PERIOD: married to high ability type
            if (chain_martm(s,tm-1) == i) then
            call random_number(rr)
            if(rr<lambda(tm-1,2,2,i-9,chain_prodmm(s,tm-1))) then  ! getting divorced
            chain_martm(s,tm) = 0
            
            ! productivity
            if (chain_agg(s,tm) == 1) then    ! boom  
            call random_number(rr)
            ZC_zc = cumsum(pic_boom(chain_prodc(s+SS,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s+SS,tm)=j 
            end do
            else 							  ! recession
            call random_number(rr)
            ZC_zc = cumsum(pic_rec(chain_prodc(s+SS,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s+SS,tm)=j 
            end do
            end if
            
            ! "translating" into male state variable
            do m=1,nz
            do j= 1,nzz
            do k= 1,nz
            
            if(chain_prodc(s+SS,tm) == (m-1)*nzz*nz+(j-1)*nz+k) then
            chain_prodm(s,tm) = j
            end if
            
            end do 
            end do
            end do      
            
            else        ! staying married
            
            ! productivity
            if (chain_agg(s,tm) == 1) then    ! boom   
            call random_number(rr)
            ZC_zc = cumsum(pic_boom(chain_prodc(s+SS,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s+SS,tm)=j 
            end do
            else 							  ! recession
            call random_number(rr)
            ZC_zc = cumsum(pic_rec(chain_prodc(s+SS,tm-1),:),nzz*nzz) 
            do j=nzz*nzz,1,-1
                if(rr<ZC_zc(j)) chain_prodc(s+SS,tm)=j 
            end do
            end if
   
            ! "translating" into male state variable
            do m=1,nz
            do j= 1,nzz
            do k= 1,nz
            
            if(chain_prodc(s+SS,tm) == (m-1)*nzz*nz+(j-1)*nz+k) then
            chain_prodmm(s,tm) = j
            chain_martm(s,tm) = (k +(m-1)*nz) + 9
            end if
            
            end do 
            end do
            end do 
            
            
            end if
            end if
            end do
            
            end if      ! end if statement whether single or married last period
            
            end do
                end do
        
        
        ! translate all back into chain_prodf and chain_prodm, for simulation purpose.
        ! in simulation part, productivity chain is defined over gender (not marital status)
        
        do s = 1, SS
                do tm = 1,TT
                
        if(chain_martm(s,tm) > 0) then
			chain_prodm(s,tm) = chain_prodmm(s,tm)
        end if    
        
        if(chain_martf(s,tm) > 0) then
			chain_prodf(s,tm) = chain_prodmf(s,tm)
        end if   
                end do
        end do
        

        ! export paths
        
        open(unit=1,file='results/chain_prodf.txt', &
        status='replace',action='write')
        open(unit=2,file='results/chain_prodm.txt', &
        status='replace',action='write')
        open(unit=3,file='results/chain_martf.txt', &
        status='replace',action='write')
        open(unit=4,file='results/chain_martm.txt', &
        status='replace',action='write')
        open(unit=5,file='results/chain_retf.txt', &
        status='replace',action='write')
        open(unit=6,file='results/chain_retm.txt', &
        status='replace',action='write')
        open(unit=7,file='results/chain_prodmf.txt', &
        status='replace',action='write')
        open(unit=8,file='results/chain_prodmm.txt', &
        status='replace',action='write')
        open(unit=9,file='results/chain_prodc.txt', &
        status='replace',action='write')
        open(unit=10,file='results/chain_agg.txt', &
        status='replace',action='write')
        
        write(1,*) chain_prodf
        write(2,*) chain_prodm
        write(3,*) chain_martf
        write(4,*) chain_martm
        write(5,*) chain_retf
        write(6,*) chain_retm
        write(7,*) chain_prodmf
        write(8,*) chain_prodmm
        write(9,*) chain_prodc
        write(10,*) chain_agg
        
        do i=1,10
        close(unit=i)
        end do
      
      write(*,*) 'shock path done'
        
end subroutine

subroutine simulate

        ! initial positions on asset grid 
        ! Note: asset is cash-on-hand, need to define lower bound in initial period
        ! lower bound is lowest possible income, dep. on family type
        
        ! couples
        do i=1,2*SS
        
        rr = rand_normal(19d0, 37d0)
        rr = max(yfc*thetafc(1)*efffc(1)*coupgrid_all_rec(1,1)+ymc*thetamc(1)*effmc(1)*coupgrid_all_rec(1,2),rr)

        grid_path_c(i,1) = cntarray(cgrid,ncc,rr)+1
        cash_path_c(i,1) = cgrid(grid_path_c(i,1))
        
        end do
        
        ! single women 
        do i=1,SS
        
		!note: mean is set lower than in data because borrowing in model not allowed, whereas data includes neg. values 
		! otherwise: do not match wealth when young: with this specification, resulting nw at 30: ca. 7, relative to 6 in data 
        rr = rand_normal(0.05d0, 16d0) 
		rr = max(yfs*thetafs(1)*efffs(1)*afs_rec(1),rr)

       ! decomposition 6
       ! rr = rand_normal(13d0, 22d0)
       ! rr = max(yfs*thetafs(1)*efffs(1)*afs_rec(1),rr)

        grid_path_f(i,1) = cntarray(cgrid,ncc,rr)+1
        cash_path_sf(i,1) = cgrid(grid_path_f(i,1))
        
        end do
        
        ! single men
        do i=1,SS
        
         
        rr = rand_normal(13d0, 22d0)
        rr = max(yms*thetams(1)*effms(1)*ams_rec(1),rr)

        grid_path_m(i,1) = cntarray(cgrid,ncc,rr)+1
        cash_path_sm(i,1)= cgrid(grid_path_m(i,1))
        
        end do
    
        
        !**********************************************************
        !************************ WOMEN ***************************
        !**********************************************************

        !!!!!!!!!!!!!!!!  low ability type  !!!!!!!!!!!!!!!!!!!!
        
        do s=1,AB_f
            do tm= 1,TT-1
            
            if(chain_martf(s,tm) == 0) then         ! CURRENT PERIOD: single 
            
            ! INTERTEMPORAL DECISIONS
            cprime_path_sf(s,tm) = cprimepol_f(grid_path_f(s,tm),chain_prodf(s,tm),1,tm)
           
            if(chain_martf(s,tm+1)==0) then         ! NEXT PERIOD: single                                 
            cash_path_sf(s,tm+1)  = yfs*efffs(tm+1)*thetafs(1)*afs_comb(chain_prodf(s,tm+1),chain_agg(s,tm+1)) + &
									((1d0-sharepol_f(grid_path_f(s,tm),chain_prodf(s,tm),1,tm))*(1d0+r) + &
									sharepol_f(grid_path_f(s,tm),chain_prodf(s,tm),1,tm)*ard_comb(chain_retf(s,tm+1),chain_agg(s,tm+1)))*cprime_path_sf(s,tm)
               
            ! where to locate on grid?  
            gridl_path_f(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_sf(s,tm+1)),1)   
            gridh_path_f(s,tm+1) = min(gridl_path_f(s,tm+1)+1,ncc)
            
            if(gridl_path_f(s,tm+1) == gridh_path_f(s,tm+1)) then
            grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else   
            
            diff = (cash_path_sf(s,tm+1)- cgrid(gridl_path_f(s,tm+1)))/(cgrid(gridh_path_f(s,tm+1))-cgrid(gridl_path_f(s,tm+1)))
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else
                grid_path_f(s,tm+1) = gridl_path_f(s,tm+1)
            end if
            end if
            
            else                                      ! NEXT PERIOD: married
            
            if(chain_martf(s,tm+1)>0 .AND. chain_martf(s,tm+1)<10 ) then   ! to low ability spouse
            
            asset_sp_ind = max(cntarray(cgrid,ncc,mwealth(1,1,tm+1)),1)  
            asset_sp = cgrid(asset_sp_ind)

            ! translate into individual income for both spouses
            if(chain_agg(s,tm+1) == 1) then    ! boom 
            if(chain_prodmf(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_boom(nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1),:)
            else if (chain_prodmf(s,tm+1) > 3 .AND. chain_prodmf(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_prodmf(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if 
            elseif (chain_agg(s,tm+1) == 2) then ! recession
            if(chain_prodmf(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_rec(nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1),:)
            else if (chain_prodmf(s,tm+1) > 3 .AND. chain_prodmf(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_prodmf(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if 
            end if

            cash_path_c(s,tm+1)  = yfc*efffc(tm+1)*thetafc(1)*coupstate(1)+ymc*effmc(tm+1)*thetamc(1)*coupstate(2) + &
								((1d0-sharepol_f(grid_path_f(s,tm),chain_prodf(s,tm),1,tm))*(1d0+r) + &
                                sharepol_f(grid_path_f(s,tm),chain_prodf(s,tm),1,tm)*ard_comb(chain_retf(s,tm+1),chain_agg(s,tm+1)))*&
                                cprime_path_sf(s,tm)+asset_sp
               
            ! where to locate on grid (cash_path)?  
            gridl_path_f(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_c(s,tm+1)),1)   
            gridh_path_f(s,tm+1) = min(gridl_path_f(s,tm+1)+1,ncc)
            
            if(gridl_path_f(s,tm+1) == gridh_path_f(s,tm+1)) then
            grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else 
            
            diff = (cash_path_c(s,tm+1)- cgrid(gridl_path_f(s,tm+1)))/(cgrid(gridh_path_f(s,tm+1))-cgrid(gridl_path_f(s,tm+1)))
            call random_number(draw)    
              if(draw<=diff) then
                grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else
                grid_path_f(s,tm+1) = gridl_path_f(s,tm+1)
            end if
            end if
            
            elseif (chain_martf(s,tm+1)>9) then			! to high ability spouse
            
            asset_sp_ind = max(cntarray(cgrid,ncc,mwealth(1,2,tm+1)),1)            
            asset_sp = cgrid(asset_sp_ind)
            
            ! translate into individual income for both spouses
            if (chain_agg(s,tm+1) == 1) then    ! boom 
            if(chain_prodmf(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_boom(nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1),:)
            else if (chain_prodmf(s,tm+1) > 3 .AND. chain_prodmf(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_prodmf(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if 
            elseif (chain_agg(s,tm+1) == 2) then ! recession
            if(chain_prodmf(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_rec(nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1),:)
            else if (chain_prodmf(s,tm+1) > 3 .AND. chain_prodmf(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_prodmf(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if 
            end if
            
            cash_path_c(s,tm+1)  = yfc*efffc(tm+1)*thetafc(1)*coupstate(1)+ymc*effmc(tm+1)*thetamc(2)*coupstate(2) + &
								((1d0-sharepol_f(grid_path_f(s,tm),chain_prodf(s,tm),1,tm))*(1d0+r) + &
                                sharepol_f(grid_path_f(s,tm),chain_prodf(s,tm),1,tm)*ard_comb(chain_retf(s,tm+1),chain_agg(s,tm+1)))*&
                                cprime_path_sf(s,tm)+asset_sp
               
            ! where to locate on grid (cash_path)?  
            gridl_path_f(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_c(s,tm+1)),1)   
            gridh_path_f(s,tm+1) = min(gridl_path_f(s,tm+1)+1,ncc)
            
            if(gridl_path_f(s,tm+1) == gridh_path_f(s,tm+1)) then
            grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else 
            
            diff = (cash_path_c(s,tm+1)- cgrid(gridl_path_f(s,tm+1)))/(cgrid(gridh_path_f(s,tm+1))-cgrid(gridl_path_f(s,tm+1)))
            call random_number(draw)    
              if(draw<=diff) then
                grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else
                grid_path_f(s,tm+1) = gridl_path_f(s,tm+1)
            end if
            end if
            
            end if
            
            
            end if 
            
            ! INTRATEMPORAL DECISIONS
            c_path_sf(s,tm)     = cpol_f(grid_path_f(s,tm),chain_prodf(s,tm),1,tm)
            share_path_sf(s,tm) = sharepol_f(grid_path_f(s,tm),chain_prodf(s,tm),1,tm)
            risky_path_sf(s,tm) = sharepol_f(grid_path_f(s,tm),chain_prodf(s,tm),1,tm)*cprime_path_sf(s,tm)
            safe_path_sf(s,tm)  = (1d0-sharepol_f(grid_path_f(s,tm),chain_prodf(s,tm),1,tm))*cprime_path_sf(s,tm)
        
            ! continuation values couple vs. single
            cntsing_path_sf(s,tm) = ct_sing_f(grid_path_f(s,tm),chain_prodf(s,tm),1,tm) 
            cntcoup_path_sf(s,tm) = ct_coup_f(grid_path_f(s,tm),chain_prodf(s,tm),1,tm)  
            
            ! human capital values    
            hcval_path_sf(s,tm)   = hc_val_f_low(tm,grid_path_f(s,tm),chain_prodf(s,tm))
  
			! stock market participation rate
            if(share_path_sf(s,tm)>0d0) then
            SMP_path_sf(s,tm) = 1
            else 
            SMP_path_sf(s,tm) = 0
            end if 

            else
            
            do i = 1, 9  !  CURRENT PERIOD: married to low ability type
            
            ! INTERTEMPORAL DECISIONS
            if(chain_martf(s,tm) == i) then         
            cprime_path_c(s,tm)  = cprimepol_c(tm,1,1,grid_path_f(s,tm),chain_prodf(s,tm),i)

            if(chain_martf(s,tm+1) > 0) then     ! NEXT PERIOD: married 
            
            ! translate into individual income for both spouses
            if (chain_agg(s,tm+1) == 1) then    ! boom 
            if(chain_prodmf(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_boom(nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1),:)
            else if (chain_prodmf(s,tm+1) > 3 .AND. chain_prodmf(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_prodmf(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if 
            elseif (chain_agg(s,tm+1) == 2) then ! recession
            if(chain_prodmf(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_rec(nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1),:)
            else if (chain_prodmf(s,tm+1) > 3 .AND. chain_prodmf(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_prodmf(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if 
            end if
                                 
            cash_path_c(s,tm+1)  = yfc*efffc(tm+1)*thetafc(1)*coupstate(1)+ymc*effmc(tm+1)*thetamc(1)*coupstate(2) + &
					   ((1d0-sharepol_c(tm,1,1,grid_path_f(s,tm),chain_prodf(s,tm),i))*(1d0+r) + &
                        sharepol_c(tm,1,1,grid_path_f(s,tm),chain_prodf(s,tm),i)*ard_comb(chain_retf(s,tm+1),chain_agg(s,tm+1)))*cprime_path_c(s,tm)
               
            ! where to locate on grid (cash_path)?  
            gridl_path_f(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_c(s,tm+1)),1)   
            gridh_path_f(s,tm+1) = min(gridl_path_f(s,tm+1)+1,ncc)
            
            if(gridl_path_f(s,tm+1)==gridh_path_f(s,tm+1)) then
            grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else
       
            diff = (cash_path_c(s,tm+1)- cgrid(gridl_path_f(s,tm+1)))/(cgrid(gridh_path_f(s,tm+1))-cgrid(gridl_path_f(s,tm+1)))
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else
                grid_path_f(s,tm+1) = gridl_path_f(s,tm+1)
            end if
            end if
   
            else                                  ! NEXT PERIOD: single (asset split 50-50)
            
            cash_path_sf(s,tm+1)  = yfs*efffs(tm+1)*thetafs(1)*afs_comb(chain_prodf(s,tm+1),chain_agg(s,tm+1)) + &
					((1d0-sharepol_c(tm,1,1,grid_path_f(s,tm),chain_prodf(s,tm),i))*(1d0+r) + &
                        sharepol_c(tm,1,1,grid_path_f(s,tm),chain_prodf(s,tm),i)* &
                        ard_comb(chain_retf(s,tm+1),chain_agg(s,tm+1)))*(cprime_path_c(s,tm)*0.45d0)
               
            ! where to locate on grid (cash_path)?  
            gridl_path_f(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_sf(s,tm+1)),1)   
            gridh_path_f(s,tm+1) = min(gridl_path_f(s,tm+1)+1,ncc)
            
            if(gridh_path_f(s,tm+1)==  gridl_path_f(s,tm+1)) then
            grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else
        
            diff = (cash_path_sf(s,tm+1)- cgrid(gridl_path_f(s,tm+1)))/(cgrid(gridh_path_f(s,tm+1))-cgrid(gridl_path_f(s,tm+1)))
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else
                grid_path_f(s,tm+1) = gridl_path_f(s,tm+1)
            end if
            end if 
            
            end if 
            
            ! INTRATEMPORAL DECISIONS
            c_path_c(s,tm)      = cpol_c(tm,1,1,grid_path_f(s,tm),chain_prodf(s,tm),i)
            share_path_c(s,tm)  = sharepol_c(tm,1,1,grid_path_f(s,tm),chain_prodf(s,tm),i)
            risky_path_c(s,tm)  = sharepol_c(tm,1,1,grid_path_f(s,tm),chain_prodf(s,tm),i)*cprime_path_c(s,tm)
            safe_path_c(s,tm)   = (1d0-sharepol_c(tm,1,1,grid_path_f(s,tm),chain_prodf(s,tm),i))*cprime_path_c(s,tm)
            
            if(share_path_c(s,tm)>0d0) then
            SMP_path_c(s,tm) = 1
            else 
            SMP_path_c(s,tm) = 0
            end if 
            
            end if 
            
            end do  ! end do loop of partner's type

            
            do i = 10, 18  !  CURRENT PERIOD: married to high ability type
            
            ! INTERTEMPORAL DECISIONS
            if(chain_martf(s,tm) == i) then
            cprime_path_c(s,tm) = cprimepol_c(tm,1,2,grid_path_f(s,tm),chain_prodf(s,tm),i-9)
            
            
            if(chain_martf(s,tm+1) > 0) then  ! NEXT PERIOD: married 
            
            ! translate into individual income for both spouses
            if (chain_agg(s,tm+1) == 1) then    ! boom 
            if(chain_prodmf(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_boom(nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1),:)
            else if (chain_prodmf(s,tm+1) > 3 .AND. chain_prodmf(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_prodmf(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if
            elseif (chain_agg(s,tm+1) == 2) then ! recession
             if(chain_prodmf(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_rec(nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1),:)
            else if (chain_prodmf(s,tm+1) > 3 .AND. chain_prodmf(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_prodmf(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if
            end if
                              
            cash_path_c(s,tm+1)  = yfc*efffc(tm+1)*thetafc(1)*coupstate(1)+ymc*effmc(tm+1)*thetamc(2)*coupstate(2) + &
						((1d0-sharepol_c(tm,1,2,grid_path_f(s,tm),chain_prodf(s,tm),i-9))*(1d0+r) + &
                        sharepol_c(tm,1,2,grid_path_f(s,tm),chain_prodf(s,tm),i-9)*ard_comb(chain_retf(s,tm+1),chain_agg(s,tm+1)))*cprime_path_c(s,tm)
               
            ! where to locate on grid (cash_path)?  
            gridl_path_f(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_c(s,tm+1)),1)   
            gridh_path_f(s,tm+1) = min(gridl_path_f(s,tm+1)+1,ncc)
               
            if(gridl_path_f(s,tm+1)== gridh_path_f(s,tm+1)) then
            grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else
       
            diff = (cash_path_c(s,tm+1)- cgrid(gridl_path_f(s,tm+1)))/(cgrid(gridh_path_f(s,tm+1))-cgrid(gridl_path_f(s,tm+1)))
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else
                grid_path_f(s,tm+1) = gridl_path_f(s,tm+1)
            end if
            end if
            
            else                                ! NEXT PERIOD: single    
    
            cash_path_sf(s,tm+1)  = yfs*efffs(tm+1)*thetafs(1)*afs_comb(chain_prodf(s,tm+1),chain_agg(s,tm+1))+ &
					((1d0-sharepol_c(tm,1,2,grid_path_f(s,tm),chain_prodf(s,tm),i-9))*(1d0+r)+&
            sharepol_c(tm,1,2,grid_path_f(s,tm),chain_prodf(s,tm),i-9)*ard_comb(chain_retf(s,tm+1),chain_agg(s,tm+1)))*(cprime_path_c(s,tm)*0.45d0)
               
            ! where to locate on grid (cash_path)?  
            gridl_path_f(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_sf(s,tm+1)),1)   
            gridh_path_f(s,tm+1) = min(gridl_path_f(s,tm+1)+1,ncc)
            
            if(gridh_path_f(s,tm+1)==gridl_path_f(s,tm+1)) then
            grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else
          
            diff = (cash_path_sf(s,tm+1)- cgrid(gridl_path_f(s,tm+1)))/(cgrid(gridh_path_f(s,tm+1))-cgrid(gridl_path_f(s,tm+1)))
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else
                grid_path_f(s,tm+1) = gridl_path_f(s,tm+1)
            end if
            end if
            
            end if 
            
            ! INTRATEMPORAL DECISIONS
            c_path_c(s,tm)      = cpol_c(tm,1,2,grid_path_f(s,tm),chain_prodf(s,tm),i-9)
            share_path_c(s,tm)  = sharepol_c(tm,1,2,grid_path_f(s,tm),chain_prodf(s,tm),i-9)
            risky_path_c(s,tm)  = sharepol_c(tm,1,2,grid_path_f(s,tm),chain_prodf(s,tm),i-9)*cprime_path_c(s,tm)
            safe_path_c(s,tm)   = (1d0-sharepol_c(tm,1,2,grid_path_f(s,tm),chain_prodf(s,tm),i-9))*cprime_path_c(s,tm)
            
            if(share_path_c(s,tm)>0d0) then
            SMP_path_c(s,tm) = 1
            else 
            SMP_path_c(s,tm) = 0
            end if 

            end if 
            end do
            
           
            end if 
            
           
            end do
        end do
        
          
        !!!!!!!!!!!!!!!!!!!!  women, high ability type  !!!!!!!!!!!!!!!!!!!!

        do s=AB_f+1,SS
            do tm= 1,TT-1
            
            if(chain_martf(s,tm) == 0) then !  CURRENT PERIOD: single
            
            ! INTERTEMPORAL DECISIONS
            cprime_path_sf(s,tm) = cprimepol_f(grid_path_f(s,tm),chain_prodf(s,tm),2,tm)
            
            if(chain_martf(s,tm+1) == 0) then     !  NEXT PERIOD: single                                
            cash_path_sf(s,tm+1)  = yfs*efffs(tm+1)*thetafs(2)*afs_comb(chain_prodf(s,tm+1),chain_agg(s,tm+1)) + &
								((1d0-sharepol_f(grid_path_f(s,tm),chain_prodf(s,tm),2,tm))*(1d0+r) + &
                                sharepol_f(grid_path_f(s,tm),chain_prodf(s,tm),2,tm)*ard_comb(chain_retf(s,tm+1),chain_agg(s,tm+1)))*cprime_path_sf(s,tm)
               
            ! where to locate on grid (cash_path)?  
            gridl_path_f(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_sf(s,tm+1)),1)   
            gridh_path_f(s,tm+1) = min(gridl_path_f(s,tm+1)+1,ncc)
            
            if(gridl_path_f(s,tm+1) == gridh_path_f(s,tm+1)) then
            grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else 
          
            diff = (cash_path_sf(s,tm+1)- cgrid(gridl_path_f(s,tm+1)))/(cgrid(gridh_path_f(s,tm+1))-cgrid(gridl_path_f(s,tm+1)))
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else
                grid_path_f(s,tm+1) = gridl_path_f(s,tm+1)
            end if
            end if
            
            else                                   ! NEXT PERIOD: married
            
            
            if(chain_martf(s,tm+1)>0 .AND. chain_martf(s,tm+1)<10 ) then  ! to low ability type
            
            asset_sp_ind = max(cntarray(cgrid,ncc,mwealth(1,1,tm+1)),1)            
            asset_sp = cgrid(asset_sp_ind)
            
            ! translate into individual income for both spouses
            if (chain_agg(s,tm+1) == 1) then    ! boom 
            if(chain_prodmf(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_boom(nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1),:)
            else if (chain_prodmf(s,tm+1) > 3 .AND. chain_prodmf(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_prodmf(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if
            elseif (chain_agg(s,tm+1) == 2) then  ! recession
             if(chain_prodmf(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_rec(nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1),:)
            else if (chain_prodmf(s,tm+1) > 3 .AND. chain_prodmf(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_prodmf(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if
            end if
            
            cash_path_c(s,tm+1)  = yfc*efffc(tm+1)*thetafc(2)*coupstate(1)+ymc*effmc(tm+1)*thetamc(1)*coupstate(2) + &
							((1d0-sharepol_f(grid_path_f(s,tm),chain_prodf(s,tm),2,tm))*(1d0+r) + &
                                sharepol_f(grid_path_f(s,tm),chain_prodf(s,tm),2,tm)*ard_comb(chain_retf(s,tm+1),chain_agg(s,tm+1)))&
                                *cprime_path_sf(s,tm)+asset_sp
               
            ! where to locate on grid (cash_path)?  
            gridl_path_f(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_c(s,tm+1)),1)   
            gridh_path_f(s,tm+1) = min(gridl_path_f(s,tm+1)+1,ncc)
            
            if(gridl_path_f(s,tm+1)==gridh_path_f(s,tm+1)) then
            grid_path_f(s,tm+1) = gridh_path_f(s,tm+1) 
            else
            diff = (cash_path_c(s,tm+1)- cgrid(gridl_path_f(s,tm+1)))/(cgrid(gridh_path_f(s,tm+1))-cgrid(gridl_path_f(s,tm+1)))
          
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else
                grid_path_f(s,tm+1) = gridl_path_f(s,tm+1)
            end if
            end if 
            
            elseif (chain_martf(s,tm+1)>9) then  ! to high ability type
            
            asset_sp_ind = max(cntarray(cgrid,ncc,mwealth(1,2,tm+1)),1)  
            asset_sp = cgrid(asset_sp_ind)
            
            ! translate into individual income for both spouses
            if (chain_agg(s,tm+1) == 1) then    ! boom 
            if(chain_prodmf(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_boom(nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1),:)
            else if (chain_prodmf(s,tm+1) > 3 .AND. chain_prodmf(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_prodmf(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if            
			elseif (chain_agg(s,tm+1) == 2) then  ! recession
			if(chain_prodmf(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_rec(nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1),:)
            else if (chain_prodmf(s,tm+1) > 3 .AND. chain_prodmf(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_prodmf(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if   
            end if

            cash_path_c(s,tm+1)  = yfc*efffc(tm+1)*thetafc(2)*coupstate(1)+ymc*effmc(tm+1)*thetamc(2)*coupstate(2) + &
								((1d0-sharepol_f(grid_path_f(s,tm),chain_prodf(s,tm),2,tm))*(1d0+r) + &
                                sharepol_f(grid_path_f(s,tm),chain_prodf(s,tm),2,tm)*ard_comb(chain_retf(s,tm+1),chain_agg(s,tm+1)))&
                                *(cprime_path_sf(s,tm))+asset_sp
               
            ! where to locate on grid (cash_path)?  
            gridl_path_f(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_c(s,tm+1)),1)   
            gridh_path_f(s,tm+1) = min(gridl_path_f(s,tm+1)+1,ncc)
            
            if(gridl_path_f(s,tm+1)==gridh_path_f(s,tm+1)) then
            grid_path_f(s,tm+1) = gridh_path_f(s,tm+1) 
            else
            diff = (cash_path_c(s,tm+1)- cgrid(gridl_path_f(s,tm+1)))/(cgrid(gridh_path_f(s,tm+1))-cgrid(gridl_path_f(s,tm+1)))
          
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else
                grid_path_f(s,tm+1) = gridl_path_f(s,tm+1)
            end if
            end if 
            
            end if
            
            end if
            
            ! INTRATEMPORAL DECISIONS
            c_path_sf(s,tm)     = cpol_f(grid_path_f(s,tm),chain_prodf(s,tm),2,tm)
            share_path_sf(s,tm) = sharepol_f(grid_path_f(s,tm),chain_prodf(s,tm),2,tm)
            risky_path_sf(s,tm) = sharepol_f(grid_path_f(s,tm),chain_prodf(s,tm),2,tm)*cprime_path_sf(s,tm)
            safe_path_sf(s,tm)  = (1d0-sharepol_f(grid_path_f(s,tm),chain_prodf(s,tm),2,tm))*cprime_path_sf(s,tm)
            
            ! continuation values couple vs. single
            cntsing_path_sf(s,tm) = ct_sing_f(grid_path_f(s,tm),chain_prodf(s,tm),2,tm) 
            cntcoup_path_sf(s,tm) = ct_coup_f(grid_path_f(s,tm),chain_prodf(s,tm),2,tm)
        
            ! human capital values   
            hcval_path_sf(s,tm)   = hc_val_f_high(tm,grid_path_f(s,tm),chain_prodf(s,tm))

			! stock market participation rate
            if(share_path_sf(s,tm)>0d0) then
            SMP_path_sf(s,tm) = 1
            else 
            SMP_path_sf(s,tm) = 0
            end if      

        
            else
            

            do i =1,9  !  CURRENT PERIOD: married to low ability type
            
            ! INTERTEMPORAL DECISIONS
            if(chain_martf(s,tm) == i) then
            cprime_path_c(s,tm) = cprimepol_c(tm,2,1,grid_path_f(s,tm),chain_prodf(s,tm),i)
                 
            if(chain_martf(s,tm+1) >0) then             ! NEXT PERIOD: married     
             
            ! translate into individual income for both spouses
            if (chain_agg(s,tm+1) == 1) then    ! boom 
            if(chain_prodmf(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_boom(nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1),:)
            else if (chain_prodmf(s,tm+1) > 3 .AND. chain_prodmf(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_prodmf(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if 
            elseif (chain_agg(s,tm+1) == 2) then ! recession
            if(chain_prodmf(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_rec(nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1),:)
            else if (chain_prodmf(s,tm+1) > 3 .AND. chain_prodmf(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_prodmf(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_martf(s,tm+1)-1)+chain_prodmf(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if 
            end if  
                                          
            cash_path_c(s,tm+1)  =  yfc*efffc(tm+1)*thetafc(2)*coupstate(1)+ymc*effmc(tm+1)*thetamc(1)*coupstate(2) + &
						((1d0-sharepol_c(tm,2,1,grid_path_f(s,tm),chain_prodf(s,tm),i))*(1d0+r)+&
                        sharepol_c(tm,2,1,grid_path_f(s,tm),chain_prodf(s,tm),i)*ard_comb(chain_retf(s,tm+1),chain_agg(s,tm+1)))*cprime_path_c(s,tm)

            ! where to locate on grid (cash_path)?  
            gridl_path_f(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_c(s,tm+1)),1)   
            gridh_path_f(s,tm+1) = min(gridl_path_f(s,tm+1)+1,ncc)
            
            if(gridl_path_f(s,tm+1)==  gridh_path_f(s,tm+1)) then
            grid_path_f(s,tm+1) = gridh_path_f(s,tm+1) 
            else
            diff = (cash_path_c(s,tm+1)- cgrid(gridl_path_f(s,tm+1)))/(cgrid(gridh_path_f(s,tm+1))-cgrid(gridl_path_f(s,tm+1)))
        
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else
                grid_path_f(s,tm+1) = gridl_path_f(s,tm+1)
            end if
            end if
            
            else                                          ! NEXT PERIOD: single (split assets 50/50)
   
            cash_path_sf(s,tm+1)  = yfs*efffs(tm+1)*thetafs(2)*afs_comb(chain_prodf(s,tm+1),chain_agg(s,tm+1)) + &
					((1d0-sharepol_c(tm,2,1,grid_path_f(s,tm),chain_prodf(s,tm),i))*(1d0+r)+&
				sharepol_c(tm,2,1,grid_path_f(s,tm),chain_prodf(s,tm),i)*ard_comb(chain_retf(s,tm+1),chain_agg(s,tm+1)))*(cprime_path_c(s,tm)*0.45d0)
               
            ! where to locate on grid (cash_path)?  
            gridl_path_f(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_sf(s,tm+1)),1)   
            gridh_path_f(s,tm+1) = min(gridl_path_f(s,tm+1)+1,ncc)
            
            if(gridl_path_f(s,tm+1)==  gridh_path_f(s,tm+1)) then
            grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)  
            else
            diff = (cash_path_sf(s,tm+1)- cgrid(gridl_path_f(s,tm+1)))/(cgrid(gridh_path_f(s,tm+1))-cgrid(gridl_path_f(s,tm+1)))
         
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else
                grid_path_f(s,tm+1) = gridl_path_f(s,tm+1)
            end if
            end if
            
            end if 
            
            ! INTRATEMPORAL DECISIONS
            c_path_c(s,tm)      = cpol_c(tm,2,1,grid_path_f(s,tm),chain_prodf(s,tm),i)
            share_path_c(s,tm)  = sharepol_c(tm,2,1,grid_path_f(s,tm),chain_prodf(s,tm),i)
            risky_path_c(s,tm)  = sharepol_c(tm,2,1,grid_path_f(s,tm),chain_prodf(s,tm),i)*cprime_path_c(s,tm)
            safe_path_c(s,tm)   = (1d0-sharepol_c(tm,2,1,grid_path_f(s,tm),chain_prodf(s,tm),i))*cprime_path_c(s,tm)
            
            if(share_path_c(s,tm)>0d0) then
            SMP_path_c(s,tm) = 1
            else 
            SMP_path_c(s,tm) = 0
            end if 
            
            
            end if 
            
            end do
            
            do i = 10, 18 !  CURRENT PERIOD: married to high ability type
            
            ! INTERTEMPORAL DECISIONS
            if(chain_martf(s,tm) == i) then
            cprime_path_c(s,tm) = cprimepol_c(tm,2,2,grid_path_f(s,tm),chain_prodf(s,tm),i-9)

              
            if(chain_martf(s,tm+1)>0) then                  !  NEXT PERIOD: married    
            
            ! translate into individual income for both spouses
            if (chain_agg(s,tm+1) == 1) then    ! boom 
            if(chain_prodmf(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_boom(nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1),:)
            else if (chain_prodmf(s,tm+1) > 3 .AND. chain_prodmf(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_prodmf(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if 
            elseif (chain_agg(s,tm+1) == 2) then ! recession 
            if(chain_prodmf(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_rec(nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1),:)
            else if (chain_prodmf(s,tm+1) > 3 .AND. chain_prodmf(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_prodmf(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_martf(s,tm+1)-10)+chain_prodmf(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if 
            end if
                                 
            cash_path_c(s,tm+1)  = yfc*efffc(tm+1)*thetafc(2)*coupstate(1)+ymc*effmc(tm+1)*thetamc(2)*coupstate(2) + &
					((1d0-sharepol_c(tm,2,2,grid_path_f(s,tm),chain_prodf(s,tm),i-9))*(1d0+r)+&
                    sharepol_c(tm,2,2,grid_path_f(s,tm),chain_prodf(s,tm),i-9)*ard_comb(chain_retf(s,tm+1),chain_agg(s,tm+1)))*cprime_path_c(s,tm)
               
            ! where to locate on grid (cash_path)?   
            gridl_path_f(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_c(s,tm+1)),1)   
            gridh_path_f(s,tm+1) = min(gridl_path_f(s,tm+1)+1,ncc)
            
            if( gridl_path_f(s,tm+1)== gridh_path_f(s,tm+1)) then
            grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else
            diff = (cash_path_c(s,tm+1)- cgrid(gridl_path_f(s,tm+1)))/(cgrid(gridh_path_f(s,tm+1))-cgrid(gridl_path_f(s,tm+1)))
       
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else
                grid_path_f(s,tm+1) = gridl_path_f(s,tm+1)
            end if
            end if
            
            else                      ! NEXT PERIOD: single (split assets 50/50)
            
            cash_path_sf(s,tm+1) = yfs*efffs(tm+1)*thetafs(2)*afs_comb(chain_prodf(s,tm+1),chain_agg(s,tm+1))+ &
						((1d0-sharepol_c(tm,2,2,grid_path_f(s,tm),chain_prodf(s,tm),i-9))*(1d0+r)+&
            sharepol_c(tm,2,2,grid_path_f(s,tm),chain_prodf(s,tm),i-9)*ard_comb(chain_retf(s,tm+1),chain_agg(s,tm+1)))*(cprime_path_c(s,tm)*0.45d0)
               
            ! where to locate on grid (cash_path)?  
            gridl_path_f(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_sf(s,tm+1)),1)   
            gridh_path_f(s,tm+1) = min(gridl_path_f(s,tm+1)+1,ncc)
               
            if(gridl_path_f(s,tm+1)==gridh_path_f(s,tm+1)) then 
            grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else  
            diff = (cash_path_sf(s,tm+1)- cgrid(gridl_path_f(s,tm+1)))/(cgrid(gridh_path_f(s,tm+1))-cgrid(gridl_path_f(s,tm+1)))
    
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_f(s,tm+1) = gridh_path_f(s,tm+1)
            else
                grid_path_f(s,tm+1) = gridl_path_f(s,tm+1)
            end if
            end if
            
            end if
            
            ! INTRATEMPORAL DECISIONS
            c_path_c(s,tm)      = cpol_c(tm,2,2,grid_path_f(s,tm),chain_prodf(s,tm),i-9)
            share_path_c(s,tm)  = sharepol_c(tm,2,2,grid_path_f(s,tm),chain_prodf(s,tm),i-9)
            risky_path_c(s,tm)  = sharepol_c(tm,2,2,grid_path_f(s,tm),chain_prodf(s,tm),i-9)*cprime_path_c(s,tm)
            safe_path_c(s,tm)   = (1d0-sharepol_c(tm,2,2,grid_path_f(s,tm),chain_prodf(s,tm),i-9))*cprime_path_c(s,tm)
            
            if(share_path_c(s,tm)>0d0) then
            SMP_path_c(s,tm) = 1
            else 
            SMP_path_c(s,tm) = 0
            end if 

            
            end if 
            end do
        
              end if
            end do
          end do
          
          
          
          !**********************************************************
          !************************** MEN ***************************
          !**********************************************************
           
          !!!!!!!!!!!!!!!!!!!!  low ability type  !!!!!!!!!!!!!!!!!!!
        
            do s=1,AB_m
            do tm= 1,TT-1
            
            if(chain_martm(s,tm) == 0) then !  CURRENT PERIOD: single
            
            ! INTERTEMPORAL DECISIONS
            cprime_path_sm(s,tm) = cprimepol_m(grid_path_m(s,tm),chain_prodm(s,tm),1,tm)

            if(chain_martm(s,tm+1)==0) then   !  NEXT PERIOD: single                                
            cash_path_sm(s,tm+1)  = yms*effms(tm+1)*thetams(1)*ams_comb(chain_prodm(s,tm+1),chain_agg(s,tm+1)) + &
				((1d0-sharepol_m(grid_path_m(s,tm),chain_prodm(s,tm),1,tm))*(1d0+r) + &
                                sharepol_m(grid_path_m(s,tm),chain_prodm(s,tm),1,tm)*ard_comb(chain_retm(s,tm+1),chain_agg(s,tm+1)))*cprime_path_sm(s,tm)
               
            ! where to locate on grid (cash_path)?  
            gridl_path_m(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_sm(s,tm+1)),1)   
            gridh_path_m(s,tm+1) = min(gridl_path_m(s,tm+1)+1,ncc)
               
            if(gridl_path_m(s,tm+1)==gridh_path_m(s,tm+1)) then
            grid_path_m(s,tm+1) = gridh_path_m(s,tm+1) 
            else
            diff = (cash_path_sm(s,tm+1)- cgrid(gridl_path_m(s,tm+1)))/(cgrid(gridh_path_m(s,tm+1))-cgrid(gridl_path_m(s,tm+1)))
        
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_m(s,tm+1) = gridh_path_m(s,tm+1)
            else
                grid_path_m(s,tm+1) = gridl_path_m(s,tm+1)
            end if
            end if 
            
            else                               !  NEXT PERIOD: married       
            
            if(chain_martm(s,tm+1)>0 .AND. chain_martm(s,tm+1)<10 ) then   ! to low ability type
            
            asset_sp_ind = max(cntarray(cgrid,ncc,mwealth(2,1,tm+1)),1)          
            asset_sp = cgrid(asset_sp_ind)
            
            ! translate into individual income within couple
            if (chain_agg(s,tm+1) == 1) then    ! boom 
            if(chain_martm(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_boom(nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1),:)
            else if (chain_martm(s,tm+1) > 3 .AND. chain_martm(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_martm(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if 
            elseif (chain_agg(s,tm+1) == 2) then ! recession
            if(chain_martm(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_rec(nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1),:)
            else if (chain_martm(s,tm+1) > 3 .AND. chain_martm(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_martm(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if 
            end if 
            
            cash_path_c(s+SS,tm+1)  = yfc*efffc(tm+1)*thetafc(1)*coupstate(1)+ ymc*effmc(tm+1)*thetamc(1)*coupstate(2) + &
						((1d0-sharepol_m(grid_path_m(s,tm),chain_prodm(s,tm),1,tm))*(1d0+r) + &
                                sharepol_m(grid_path_m(s,tm),chain_prodm(s,tm),1,tm)*ard_comb(chain_retm(s,tm+1),chain_agg(s,tm+1)))&
                                *cprime_path_sm(s,tm)+asset_sp
               
            ! where to locate on grid (cash_path)?  
            gridl_path_m(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_c(s+SS,tm+1)),1)   
            gridh_path_m(s,tm+1) = min(gridl_path_m(s,tm+1)+1,ncc)
            
            if( gridl_path_m(s,tm+1)== gridh_path_m(s,tm+1)) then
            grid_path_m(s,tm+1) = gridh_path_m(s,tm+1) 
            else  
            diff = (cash_path_c(s+SS,tm+1)- cgrid(gridl_path_m(s,tm+1)))/(cgrid(gridh_path_m(s,tm+1))-cgrid(gridl_path_m(s,tm+1)))
          
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_m(s,tm+1) = gridh_path_m(s,tm+1)
            else
                grid_path_m(s,tm+1) = gridl_path_m(s,tm+1)
            end if
            end if
            
            elseif (chain_martm(s,tm+1)>9) then  ! to high ability type
            
            asset_sp_ind = max(cntarray(cgrid,ncc,mwealth(2,2,tm+1)),1)           
            asset_sp = cgrid(asset_sp_ind)
            
            ! translate into individual income within couple
            if (chain_agg(s,tm+1) == 1) then    ! boom 
            if(chain_martm(s,tm+1)-9 <= 3) then
              coupstate(:) = coupgrid_all_boom(nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9,:)
            else if (chain_martm(s,tm+1)-9 > 3 .AND. chain_martm(s,tm+1)-9 < 7) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9-nz)+nz*nzz,:)
            else if  (chain_martm(s,tm+1)-9 > 6) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9-2*nz)+2*nz*nzz,:)
            end if 
            elseif (chain_agg(s,tm+1) == 2) then ! recession
            if(chain_martm(s,tm+1)-9 <= 3) then
              coupstate(:) = coupgrid_all_rec(nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9,:)
            else if (chain_martm(s,tm+1)-9 > 3 .AND. chain_martm(s,tm+1)-9 < 7) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9-nz)+nz*nzz,:)
            else if  (chain_martm(s,tm+1)-9 > 6) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9-2*nz)+2*nz*nzz,:)
            end if 
            end if
            
            cash_path_c(s+SS,tm+1)  = yfc*efffc(tm+1)*thetafc(2)*coupstate(1)+ ymc*effmc(tm+1)*thetamc(1)*coupstate(2) + &
									((1d0-sharepol_m(grid_path_m(s,tm),chain_prodm(s,tm),1,tm))*(1d0+r) + &
                                sharepol_m(grid_path_m(s,tm),chain_prodm(s,tm),1,tm)*ard_comb(chain_retm(s,tm+1),chain_agg(s,tm+1)))&
                                *(cprime_path_sm(s,tm))+asset_sp
               
            ! where to locate on grid (cash_path)?  
            gridl_path_m(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_c(s+SS,tm+1)),1)   
            gridh_path_m(s,tm+1) = min(gridl_path_m(s,tm+1)+1,ncc)
            
            if( gridl_path_m(s,tm+1)== gridh_path_m(s,tm+1)) then
            grid_path_m(s,tm+1) = gridh_path_m(s,tm+1) 
            else  
            diff = (cash_path_c(s+SS,tm+1)- cgrid(gridl_path_m(s,tm+1)))/(cgrid(gridh_path_m(s,tm+1))-cgrid(gridl_path_m(s,tm+1)))
          
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_m(s,tm+1) = gridh_path_m(s,tm+1)
            else
                grid_path_m(s,tm+1) = gridl_path_m(s,tm+1)
            end if
            end if
            
            end if
            
            end if 
            
            ! INTRATEMPORAL DECISIONS
            c_path_sm(s,tm)     = cpol_m(grid_path_m(s,tm),chain_prodm(s,tm),1,tm)
            share_path_sm(s,tm) = sharepol_m(grid_path_m(s,tm),chain_prodm(s,tm),1,tm)
            risky_path_sm(s,tm) = sharepol_m(grid_path_m(s,tm),chain_prodm(s,tm),1,tm)*cprime_path_sm(s,tm)
            safe_path_sm(s,tm)  = (1d0-sharepol_m(grid_path_m(s,tm),chain_prodm(s,tm),1,tm))*cprime_path_sm(s,tm)
            
            ! continuation values couple vs. single
            cntsing_path_sm(s,tm) = ct_sing_m(grid_path_m(s,tm),chain_prodm(s,tm),1,tm) 
            cntcoup_path_sm(s,tm) = ct_coup_m(grid_path_m(s,tm),chain_prodm(s,tm),1,tm)
            
            ! human capital values    
            hcval_path_sm(s,tm)   = hc_val_m_low(tm,grid_path_m(s,tm),chain_prodm(s,tm))

            ! stock market participation rate
            if(share_path_sm(s,tm)>0d0) then
            SMP_path_sm(s,tm) = 1
            else 
            SMP_path_sm(s,tm) = 0
            end if 

            else
            
            
            do i = 1, 9  !  CURRENT PERIOD: married to low ability type
            
            ! INTERTEMPORAL DECISIONS
            if(chain_martm(s,tm) == i) then
            cprime_path_c(s+SS,tm) = cprimepol_c(tm,1,1,grid_path_m(s,tm),i,chain_prodm(s,tm))

            if(chain_martm(s,tm+1) > 0) then    !  NEXT PERIOD: married 
            
            ! translate into individual income within couple
            if (chain_agg(s,tm+1) == 1) then    ! boom
            if(chain_martm(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_boom(nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1),:)
            else if (chain_martm(s,tm+1) > 3 .AND. chain_martm(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_martm(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if 
            elseif (chain_agg(s,tm+1) == 2) then  ! recession
            if(chain_martm(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_rec(nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1),:)
            else if (chain_martm(s,tm+1) > 3 .AND. chain_martm(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_martm(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if 
            end if
                                              
            cash_path_c(s+SS,tm+1)  = yfc*efffc(tm+1)*thetafc(1)*coupstate(1)+ ymc*effmc(tm+1)*thetamc(1)*coupstate(2) + &
						((1d0-sharepol_c(tm,1,1,grid_path_m(s,tm),i,chain_prodm(s,tm)))*(1d0+r)+&
                        sharepol_c(tm,1,1,grid_path_m(s,tm),i,chain_prodm(s,tm))*ard_comb(chain_retm(s,tm+1),chain_agg(s,tm+1)))*cprime_path_c(s+SS,tm)
               
            ! where to locate on grid (cash_path)?  
            gridl_path_m(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_c(s+SS,tm+1)),1)   
            gridh_path_m(s,tm+1) = min(gridl_path_m(s,tm+1)+1,ncc)
            if(gridl_path_m(s,tm+1)==gridh_path_m(s,tm+1)) then
            grid_path_m(s,tm+1) = gridh_path_m(s,tm+1)
            else
            diff = (cash_path_c(s+SS,tm+1)- cgrid(gridl_path_m(s,tm+1)))/(cgrid(gridh_path_m(s,tm+1))-cgrid(gridl_path_m(s,tm+1)))
       
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_m(s,tm+1) = gridh_path_m(s,tm+1)
            else
                grid_path_m(s,tm+1) = gridl_path_m(s,tm+1)
            end if
            end if
   
            else                                ! NEXT PERIOD: single (split assets 50/50)
            
            cash_path_sm(s,tm+1)  = yms*effms(tm+1)*thetams(1)*ams_comb(chain_prodm(s,tm+1),chain_agg(s,tm+1)) + &
					((1d0-sharepol_c(tm,1,1,grid_path_m(s,tm),i,chain_prodm(s,tm)))*(1d0+r) + &
				sharepol_c(tm,1,1,grid_path_m(s,tm),i,chain_prodm(s,tm))*ard_comb(chain_retm(s,tm+1),chain_agg(s,tm+1)))*(cprime_path_c(s+SS,tm)*0.45d0)
               
            ! where to locate on grid (cash_path)?  
            gridl_path_m(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_sm(s,tm+1)),1)   
            gridh_path_m(s,tm+1) = min(gridl_path_m(s,tm+1)+1,ncc)
            
            if(gridl_path_m(s,tm+1)==gridh_path_m(s,tm+1)) then
            grid_path_m(s,tm+1) = gridh_path_m(s,tm+1)
            else   
            diff = (cash_path_sm(s,tm+1)- cgrid(gridl_path_m(s,tm+1)))/(cgrid(gridh_path_m(s,tm+1))-cgrid(gridl_path_m(s,tm+1)))
     
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_m(s,tm+1) = gridh_path_m(s,tm+1)
            else
                grid_path_m(s,tm+1) = gridl_path_m(s,tm+1)
            end if
            end if
            
            
            end if 
            
            ! INTRATEMPORAL DECISIONS
            c_path_c(s+SS,tm)      = cpol_c(tm,1,1,grid_path_m(s,tm),i,chain_prodm(s,tm))
            share_path_c(s+SS,tm)  = sharepol_c(tm,1,1,grid_path_m(s,tm),i,chain_prodm(s,tm))
            risky_path_c(s+SS,tm)  = sharepol_c(tm,1,1,grid_path_m(s,tm),i,chain_prodm(s,tm))*cprime_path_c(s+SS,tm)
            safe_path_c(s+SS,tm)   = (1d0-sharepol_c(tm,1,1,grid_path_m(s,tm),i,chain_prodm(s,tm)))*cprime_path_c(s+SS,tm)
            
            if(share_path_c(s+SS,tm)>0d0) then
            SMP_path_c(s+SS,tm) = 1
            else 
            SMP_path_c(s+SS,tm) = 0
            end if 
            
            end if 
            end do
            
            do i = 10, 18  !  CURRENT PERIOD: married to high ability type
            
            ! INTERTEMPORAL DECISIONS
            if(chain_martm(s,tm) == i) then
            cprime_path_c(s+SS,tm) = cprimepol_c(tm,2,1,grid_path_m(s,tm),i-9,chain_prodm(s,tm))
            
            if(chain_martm(s,tm+1) > 0) then  ! NEXT PERIOD: married
            
            ! translate into individual income within couple
            if (chain_agg(s,tm+1) == 1) then    ! boom
            if(chain_martm(s,tm+1)-9 <= 3) then
              coupstate(:) = coupgrid_all_boom(nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9,:)
            else if (chain_martm(s,tm+1)-9 > 3 .AND. chain_martm(s,tm+1)-9 < 7) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9-nz)+nz*nzz,:)
            else if  (chain_martm(s,tm+1)-9 > 6) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9-2*nz)+2*nz*nzz,:)
            end if 
            elseif (chain_agg(s,tm+1) == 2) then ! recession
            if(chain_martm(s,tm+1)-9 <= 3) then
              coupstate(:) = coupgrid_all_rec(nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9,:)
            else if (chain_martm(s,tm+1)-9 > 3 .AND. chain_martm(s,tm+1)-9 < 7) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9-nz)+nz*nzz,:)
            else if  (chain_martm(s,tm+1)-9 > 6) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9-2*nz)+2*nz*nzz,:)
            end if
            end if
                                             
            cash_path_c(s+SS,tm+1)  = yfc*efffc(tm+1)*thetafc(2)*coupstate(1)+ ymc*effmc(tm+1)*thetamc(1)*coupstate(2) + &
					((1d0-sharepol_c(tm,2,1,grid_path_m(s,tm),i-9,chain_prodm(s,tm)))*(1d0+r) + &
                        sharepol_c(tm,2,1,grid_path_m(s,tm),i-9,chain_prodm(s,tm))*ard_comb(chain_retm(s,tm+1),chain_agg(s,tm+1)))*cprime_path_c(s+SS,tm)
               
            ! where to locate on grid (cash_path)?  
            gridl_path_m(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_c(s+SS,tm+1)),1)   
            gridh_path_m(s,tm+1) = min(gridl_path_m(s,tm+1)+1,ncc)
            
            if(gridl_path_m(s,tm+1)==gridh_path_m(s,tm+1)) then
            grid_path_m(s,tm+1) = gridh_path_m(s,tm+1) 
            else 
            diff = (cash_path_c(s+SS,tm+1)- cgrid(gridl_path_m(s,tm+1)))/(cgrid(gridh_path_m(s,tm+1))-cgrid(gridl_path_m(s,tm+1)))
         
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_m(s,tm+1) = gridh_path_m(s,tm+1)
            else
                grid_path_m(s,tm+1) = gridl_path_m(s,tm+1)
            end if
            end if 
            
            else                              ! NEXT PERIOD: single 
    
            cash_path_sm(s,tm+1)  = yms*effms(tm+1)*thetams(1)*ams_comb(chain_prodm(s,tm+1),chain_agg(s,tm+1)) + &
					((1d0-sharepol_c(tm,2,1,grid_path_m(s,tm),i-9,chain_prodm(s,tm)))*(1d0+r) + &
				sharepol_c(tm,2,1,grid_path_m(s,tm),i-9,chain_prodm(s,tm))*ard_comb(chain_retm(s,tm+1),chain_agg(s,tm+1)))*(cprime_path_c(s+SS,tm)*0.45d0)
               
            ! where to locate on grid (cash_path)?  
            gridl_path_m(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_sm(s,tm+1)),1)   
            gridh_path_m(s,tm+1) = min(gridl_path_m(s,tm+1)+1,ncc)
             
            if(gridl_path_m(s,tm+1)==gridh_path_m(s,tm+1)) then
            grid_path_m(s,tm+1) = gridh_path_m(s,tm+1) 
            else    
            diff = (cash_path_sm(s,tm+1)- cgrid(gridl_path_m(s,tm+1)))/(cgrid(gridh_path_m(s,tm+1))-cgrid(gridl_path_m(s,tm+1)))
          
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_m(s,tm+1) = gridh_path_m(s,tm+1)
            else
                grid_path_m(s,tm+1) = gridl_path_m(s,tm+1)
            end if
            end if
            
            end if 
            
            ! INTRATEMPORAL DECISIONS
            c_path_c(s+SS,tm)      = cpol_c(tm,2,1,grid_path_m(s,tm),i-9,chain_prodm(s,tm))
            share_path_c(s+SS,tm)  = sharepol_c(tm,2,1,grid_path_m(s,tm),i-9,chain_prodm(s,tm))
            risky_path_c(s+SS,tm)  = sharepol_c(tm,2,1,grid_path_m(s,tm),i-9,chain_prodm(s,tm))*cprime_path_c(s+SS,tm)
            safe_path_c(s+SS,tm)   = (1d0-sharepol_c(tm,2,1,grid_path_m(s,tm),i-9,chain_prodm(s,tm)))*cprime_path_c(s+SS,tm)
            
            
            if(share_path_c(s+SS,tm)>0d0) then
            SMP_path_c(s+SS,tm) = 1
            else 
            SMP_path_c(s+SS,tm) = 0
            end if 
            
            
            end if 
            end do
            
           
            end if 
           
        
            end do
        end do
        

         !!!!!!!!!!!!!!!!!!!! high ability type  !!!!!!!!!!!!!!!!!!!!

        do s=AB_m+1,SS
            do tm= 1,TT-1

            
            if(chain_martm(s,tm) == 0) then   !  CURRENT PERIOD: single
            
            ! INTERTEMPORAL DECISIONS
            cprime_path_sm(s,tm) = cprimepol_m(grid_path_m(s,tm),chain_prodm(s,tm),2,tm)
            
            if(chain_martm(s,tm+1) == 0) then     !  NEXT PERIOD: single                        
            cash_path_sm(s,tm+1)  = yms*effms(tm+1)*thetams(2)*ams_comb(chain_prodm(s,tm+1),chain_agg(s,tm+1)) + &
						((1d0-sharepol_m(grid_path_m(s,tm),chain_prodm(s,tm),2,tm))*(1d0+r) + &
                                sharepol_m(grid_path_m(s,tm),chain_prodm(s,tm),2,tm)*ard_comb(chain_retm(s,tm+1),chain_agg(s,tm+1)))*cprime_path_sm(s,tm)
               
            ! where to locate on grid (cash_path)?  
            gridl_path_m(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_sm(s,tm+1)),1)   
            gridh_path_m(s,tm+1) = min(gridl_path_m(s,tm+1)+1,ncc)
            
            if(gridl_path_m(s,tm+1)==gridh_path_m(s,tm+1)) then
            grid_path_m(s,tm+1) = gridh_path_m(s,tm+1)
            else   
            diff = (cash_path_sm(s,tm+1)- cgrid(gridl_path_m(s,tm+1)))/(cgrid(gridh_path_m(s,tm+1))-cgrid(gridl_path_m(s,tm+1)))
      
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_m(s,tm+1) = gridh_path_m(s,tm+1)
            else
                grid_path_m(s,tm+1) = gridl_path_m(s,tm+1)
            end if
            end if
            
            else                                   !  NEXT PERIOD: married  
            
            if(chain_martm(s,tm+1)>0 .AND. chain_martm(s,tm+1)<10 ) then   ! to low ability type
            
            asset_sp_ind = max(cntarray(cgrid,ncc,mwealth(2,1,tm+1)),1)             
            asset_sp = cgrid(asset_sp_ind)
            
		    ! translate into individual income within couple
		    if (chain_agg(s,tm+1) == 1) then    ! boom
            if(chain_martm(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_boom(nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1),:)
            else if (chain_martm(s,tm+1) > 3 .AND. chain_martm(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_martm(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if 		
            elseif (chain_agg(s,tm+1) == 2) then ! recession
            if(chain_martm(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_rec(nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1),:)
            else if (chain_martm(s,tm+1) > 3 .AND. chain_martm(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_martm(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if 	
            end if
            
            cash_path_c(s+SS,tm+1)  = yfc*efffc(tm+1)*thetafc(1)*coupstate(1)+ ymc*effmc(tm+1)*thetamc(2)*coupstate(2) + &
				((1d0-sharepol_m(grid_path_m(s,tm),chain_prodm(s,tm),2,tm))*(1d0+r) + &
                                sharepol_m(grid_path_m(s,tm),chain_prodm(s,tm),2,tm)*ard_comb(chain_retm(s,tm+1),chain_agg(s,tm+1)))*&
                                (cprime_path_sm(s,tm))+asset_sp
               
            ! where to locate on grid (cash_path)?  
            gridl_path_m(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_c(s+SS,tm+1)),1)   
            gridh_path_m(s,tm+1) = min(gridl_path_m(s,tm+1)+1,ncc)
               
            if(gridl_path_m(s,tm+1)==gridh_path_m(s,tm+1)) then
            grid_path_m(s,tm+1) = gridh_path_m(s,tm+1) 
            else
            diff = (cash_path_c(s+SS,tm+1)- cgrid(gridl_path_m(s,tm+1)))/(cgrid(gridh_path_m(s,tm+1))-cgrid(gridl_path_m(s,tm+1)))
           
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_m(s,tm+1) = gridh_path_m(s,tm+1)
            else
                grid_path_m(s,tm+1) = gridl_path_m(s,tm+1)
            end if
            end if
            
            elseif (chain_martm(s,tm+1)>9) then   ! to high ability type
            
            asset_sp_ind = max(cntarray(cgrid,ncc,mwealth(2,2,tm+1)),1)             
            asset_sp = cgrid(asset_sp_ind)
            
            ! translate into individual income within couple
            if (chain_agg(s,tm+1) == 1) then    ! boom
            if(chain_martm(s,tm+1)-9 <= 3) then
              coupstate(:) = coupgrid_all_boom(nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9,:)
            else if (chain_martm(s,tm+1)-9 > 3 .AND. chain_martm(s,tm+1)-9 < 7) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9-nz)+nz*nzz,:)
            else if  (chain_martm(s,tm+1)-9 > 6) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9-2*nz)+2*nz*nzz,:)
            end if 
            elseif (chain_agg(s,tm+1) == 2) then ! recession
            if(chain_martm(s,tm+1)-9 <= 3) then
              coupstate(:) = coupgrid_all_rec(nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9,:)
            else if (chain_martm(s,tm+1)-9 > 3 .AND. chain_martm(s,tm+1)-9 < 7) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9-nz)+nz*nzz,:)
            else if  (chain_martm(s,tm+1)-9 > 6) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9-2*nz)+2*nz*nzz,:)
            end if 
            end if
 
             cash_path_c(s+SS,tm+1)  = yfc*efffc(tm+1)*thetafc(2)*coupstate(1)+ ymc*effmc(tm+1)*thetamc(2)*coupstate(2) + &
							((1d0-sharepol_m(grid_path_m(s,tm),chain_prodm(s,tm),2,tm))*(1d0+r) + &
                                sharepol_m(grid_path_m(s,tm),chain_prodm(s,tm),2,tm)*ard_comb(chain_retm(s,tm+1),chain_agg(s,tm+1)))*&
                                cprime_path_sm(s,tm)+asset_sp
               
            ! where to locate on grid (cash_path)?  
            gridl_path_m(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_c(s+SS,tm+1)),1)   
            gridh_path_m(s,tm+1) = min(gridl_path_m(s,tm+1)+1,ncc)
               
            if(gridl_path_m(s,tm+1)==gridh_path_m(s,tm+1)) then
            grid_path_m(s,tm+1) = gridh_path_m(s,tm+1) 
            else
            diff = (cash_path_c(s+SS,tm+1)- cgrid(gridl_path_m(s,tm+1)))/(cgrid(gridh_path_m(s,tm+1))-cgrid(gridl_path_m(s,tm+1)))
           
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_m(s,tm+1) = gridh_path_m(s,tm+1)
            else
                grid_path_m(s,tm+1) = gridl_path_m(s,tm+1)
            end if
            end if
            
            end if

            end if
            
            ! INTRATEMPORAL DECISIONS
            c_path_sm(s,tm)     = cpol_m(grid_path_m(s,tm),chain_prodm(s,tm),2,tm)
            share_path_sm(s,tm) = sharepol_m(grid_path_m(s,tm),chain_prodm(s,tm),2,tm)
            risky_path_sm(s,tm) = sharepol_m(grid_path_m(s,tm),chain_prodm(s,tm),2,tm)*cprime_path_sm(s,tm)
            safe_path_sm(s,tm)  = (1d0-sharepol_m(grid_path_m(s,tm),chain_prodm(s,tm),2,tm))*cprime_path_sm(s,tm)
            
            ! continuation values couple vs. single
            cntsing_path_sm(s,tm) = ct_sing_m(grid_path_m(s,tm),chain_prodm(s,tm),2,tm) 
            cntcoup_path_sm(s,tm) = ct_coup_m(grid_path_m(s,tm),chain_prodm(s,tm),2,tm)

            ! human capital values    
            hcval_path_sm(s,tm)   = hc_val_m_high(tm,grid_path_m(s,tm),chain_prodm(s,tm))

			! stock market participation rate 
            if(share_path_sm(s,tm)>0d0) then
            SMP_path_sm(s,tm) = 1
            else 
            SMP_path_sm(s,tm) = 0
            end if      
        
            else
            
            
            do i =1,9   !  CURRENT PERIOD: married to low ability type
            
            ! INTERTEMPORAL DECISIONS
            if(chain_martm(s,tm) == i) then
            cprime_path_c(s+SS,tm) = cprimepol_c(tm,1,2,grid_path_m(s,tm),i,chain_prodm(s,tm))
                 
            if(chain_martm(s,tm+1) >0) then                     !  NEXT PERIOD: married 
            
            ! translate into individual income within couple
            if (chain_agg(s,tm+1) == 1) then    ! boom
            if(chain_martm(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_boom(nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1),:)
            else if (chain_martm(s,tm+1) > 3 .AND. chain_martm(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_martm(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if 
            elseif (chain_agg(s,tm+1) == 2) then ! recession
            if(chain_martm(s,tm+1) <= 3) then
              coupstate(:) = coupgrid_all_rec(nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1),:)
            else if (chain_martm(s,tm+1) > 3 .AND. chain_martm(s,tm+1) < 7) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-nz)+nz*nzz,:)
            else if  (chain_martm(s,tm+1) > 6) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-2*nz)+2*nz*nzz,:)
            end if 
            end if 
              
            cash_path_c(s+SS,tm+1)  = yfc*efffc(tm+1)*thetafc(1)*coupstate(1)+ ymc*effmc(tm+1)*thetamc(2)*coupstate(2) + & 
						((1d0-sharepol_c(tm,1,2,grid_path_m(s,tm),i,chain_prodm(s,tm)))*(1d0+r) + &
                            sharepol_c(tm,1,2,grid_path_m(s,tm),i,chain_prodm(s,tm))*ard_comb(chain_retm(s,tm+1),chain_agg(s,tm+1)))*cprime_path_c(s+SS,tm)
                            
            ! where to locate on grid (cash_path)?  
            gridl_path_m(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_c(s+SS,tm+1)),1)   
            gridh_path_m(s,tm+1) = min(gridl_path_m(s,tm+1)+1,ncc)
               
            if(gridl_path_m(s,tm+1)== gridh_path_m(s,tm+1)) then
            grid_path_m(s,tm+1) = gridh_path_m(s,tm+1)
            else  
            diff = (cash_path_c(s+SS,tm+1)- cgrid(gridl_path_m(s,tm+1)))/(cgrid(gridh_path_m(s,tm+1))-cgrid(gridl_path_m(s,tm+1)))
          
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_m(s,tm+1) = gridh_path_m(s,tm+1)
            else
                grid_path_m(s,tm+1) = gridl_path_m(s,tm+1)
            end if
            end if
            
            else                                                 !  NEXT PERIOD: single     
   
            cash_path_sm(s,tm+1)  = yms*effms(tm+1)*thetams(2)*ams_comb(chain_prodm(s,tm+1),chain_agg(s,tm+1)) + &
						((1d0-sharepol_c(tm,1,2,grid_path_m(s,tm),i,chain_prodm(s,tm)))*(1d0+r) + &
              sharepol_c(tm,1,2,grid_path_m(s,tm),i,chain_prodm(s,tm))*ard_comb(chain_retm(s,tm+1),chain_agg(s,tm+1)))*(cprime_path_c(s+SS,tm)*0.45d0)
               
            ! where to locate on grid (cash_path)?  
            gridl_path_m(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_sm(s,tm+1)),1)   
            gridh_path_m(s,tm+1) = min(gridl_path_m(s,tm+1)+1,ncc)
              
            if(gridl_path_m(s,tm+1)== gridh_path_m(s,tm+1)) then
            grid_path_m(s,tm+1) = gridh_path_m(s,tm+1)
            else     
            diff = (cash_path_sm(s,tm+1)- cgrid(gridl_path_m(s,tm+1)))/(cgrid(gridh_path_m(s,tm+1))-cgrid(gridl_path_m(s,tm+1)))
           
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_m(s,tm+1) = gridh_path_m(s,tm+1)
            else
                grid_path_m(s,tm+1) = gridl_path_m(s,tm+1)
            end if
            end if
            
            end if 
            
            ! INTRATEMPORAL DECISIONS
            c_path_c(s+SS,tm)      = cpol_c(tm,1,2,grid_path_m(s,tm),i,chain_prodm(s,tm))
            share_path_c(s+SS,tm)  = sharepol_c(tm,1,2,grid_path_m(s,tm),i,chain_prodm(s,tm))
            risky_path_c(s+SS,tm)  = sharepol_c(tm,1,2,grid_path_m(s,tm),i,chain_prodm(s,tm))*cprime_path_c(s+SS,tm)
            safe_path_c(s+SS,tm)   = (1d0-sharepol_c(tm,1,2,grid_path_m(s,tm),i,chain_prodm(s,tm)))*cprime_path_c(s+SS,tm)
            
            if(share_path_c(s+SS,tm)>0d0) then
            SMP_path_c(s+SS,tm) = 1
            else 
            SMP_path_c(s+SS,tm) = 0
            end if 
            
            end if 
            
            end do
            
            
            
            do i = 10, 18  !  CURRENT PERIOD: married to high ability type
            
            ! INTERTEMPORAL DECISIONS
            if(chain_martm(s,tm) == i) then
            cprime_path_c(s+SS,tm) = cprimepol_c(tm,2,2,grid_path_m(s,tm),i-9,chain_prodm(s,tm))

            if(chain_martm(s,tm+1)>0) then                  !  NEXT PERIOD: married             
            
            ! translate into individual income within couple
            if (chain_agg(s,tm+1) == 1) then    ! boom
            if(chain_martm(s,tm+1)-9 <= 3) then
              coupstate(:) = coupgrid_all_boom(nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9,:)
            else if (chain_martm(s,tm+1)-9 > 3 .AND. chain_martm(s,tm+1)-9 < 7) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9-nz)+nz*nzz,:)
            else if  (chain_martm(s,tm+1)-9 > 6) then
              coupstate(:) = coupgrid_all_boom((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9-2*nz)+2*nz*nzz,:)
            end if 
            elseif (chain_agg(s,tm+1) == 2) then ! recession
             if(chain_martm(s,tm+1)-9 <= 3) then
              coupstate(:) = coupgrid_all_rec(nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9,:)
            else if (chain_martm(s,tm+1)-9 > 3 .AND. chain_martm(s,tm+1)-9 < 7) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9-nz)+nz*nzz,:)
            else if  (chain_martm(s,tm+1)-9 > 6) then
              coupstate(:) = coupgrid_all_rec((nz*(chain_prodmm(s,tm+1)-1)+chain_martm(s,tm+1)-9-2*nz)+2*nz*nzz,:)
            end if 
            end if
                                         
            cash_path_c(s+SS,tm+1)  =  yfc*efffc(tm+1)*thetafc(2)*coupstate(1)+ ymc*effmc(tm+1)*thetamc(2)*coupstate(2) + & 
						((1d0-sharepol_c(tm,2,2,grid_path_m(s,tm),i-9,chain_prodm(s,tm)))*(1d0+r) + &
                    sharepol_c(tm,2,2,grid_path_m(s,tm),i-9,chain_prodm(s,tm))*ard_comb(chain_retm(s,tm+1),chain_agg(s,tm+1)))*cprime_path_c(s+SS,tm)
               
            ! where to locate on grid (cash_path)?  
            gridl_path_m(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_c(s+SS,tm+1)),1)   
            gridh_path_m(s,tm+1) = min(gridl_path_m(s,tm+1)+1,ncc)
            
            if(gridl_path_m(s,tm+1)== gridh_path_m(s,tm+1)) then
            grid_path_m(s,tm+1) = gridh_path_m(s,tm+1)
            else  
               
            diff = (cash_path_c(s+SS,tm+1)- cgrid(gridl_path_m(s,tm+1)))/(cgrid(gridh_path_m(s,tm+1))-cgrid(gridl_path_m(s,tm+1)))

            call random_number(draw)    
            if(draw<=diff) then
                grid_path_m(s,tm+1) = gridh_path_m(s,tm+1)
            else
                grid_path_m(s,tm+1) = gridl_path_m(s,tm+1)
            end if
            end if
            
            else                                             !  NEXT PERIOD: single     
            
            cash_path_sm(s,tm+1)  = yms*effms(tm+1)*thetams(2)*ams_comb(chain_prodm(s,tm+1),chain_agg(s,tm+1)) + &
					((1d0-sharepol_c(tm,2,2,grid_path_m(s,tm),i-9,chain_prodm(s,tm)))*(1d0+r) + &
            sharepol_c(tm,2,2,grid_path_m(s,tm),i-9,chain_prodm(s,tm))*ard_comb(chain_retm(s,tm+1),chain_agg(s,tm+1)))*(cprime_path_c(s+SS,tm)*0.45d0)
               
            ! where to locate on grid (cash_path)?  
            gridl_path_m(s,tm+1) = max(cntarray(cgrid,ncc,cash_path_sm(s,tm+1)),1)   
            gridh_path_m(s,tm+1) = min(gridl_path_m(s,tm+1)+1,ncc)
               
            if(gridl_path_m(s,tm+1)== gridh_path_m(s,tm+1)) then
            grid_path_m(s,tm+1) = gridh_path_m(s,tm+1)
            else     
               
            diff = (cash_path_sm(s,tm+1)- cgrid(gridl_path_m(s,tm+1)))/(cgrid(gridh_path_m(s,tm+1))-cgrid(gridl_path_m(s,tm+1)))
         
            call random_number(draw)    
            if(draw<=diff) then
                grid_path_m(s,tm+1) = gridh_path_m(s,tm+1)
            else
                grid_path_m(s,tm+1) = gridl_path_m(s,tm+1)
            end if
            end if
            
            end if
            
            ! INTRATEMPORAL DECISIONS   
            c_path_c(s+SS,tm)      = cpol_c(tm,2,2,grid_path_m(s,tm),i-9,chain_prodm(s,tm))
            share_path_c(s+SS,tm)  = sharepol_c(tm,2,2,grid_path_m(s,tm),i-9,chain_prodm(s,tm))
            risky_path_c(s+SS,tm)  = sharepol_c(tm,2,2,grid_path_m(s,tm),i-9,chain_prodm(s,tm))*cprime_path_c(s+SS,tm)
            safe_path_c(s+SS,tm)   = (1d0-sharepol_c(tm,2,2,grid_path_m(s,tm),i-9,chain_prodm(s,tm)))*cprime_path_c(s+SS,tm)
            
            if(share_path_c(s+SS,tm)>0d0) then
            SMP_path_c(s+SS,tm) = 1
            else 
            SMP_path_c(s+SS,tm) = 0
            end if 
            
            end if 
            end do
        
              end if
            end do
          end do
       
       
       write(*,*) 'simulation done'
          
     
end subroutine

subroutine cohort_vars
          
          ! create cohort variables on simulated data
          
          ! first, take care of values that are not allocated:
          ! that is, when single for couple and vice versa
          
          ! WOMEN
          do s=1,SS
            do tm= 1,TT-1

            inf = 900000000000d0
                      
            if(chain_martf(s,tm) > 0) then
            
            cash_path_sf(s,tm)    = -inf
            c_path_sf(s,tm)       = -inf
            SMP_path_sf(s,tm)     = -inf
            risky_path_sf(s,tm)   = -inf
            safe_path_sf(s,tm)    = -inf
            share_path_sf(s,tm)   = -inf
            hcval_path_sf(s,tm)   = -inf
            else            
            
            cash_path_c(s,tm)    = -inf
            c_path_c(s,tm)       = -inf
            SMP_path_c(s,tm)     = -inf
            risky_path_c(s,tm)   = -inf
            safe_path_c(s,tm)    = -inf
            share_path_c(s,tm)   = -inf
            
            end if
          
             end do 
            end do
          
            ! MEN
            do s=1,SS
             do tm= 1,TT-1
             
            inf = 900000000000d0
            
            if(chain_martm(s,tm) > 0) then
            
            cash_path_sm(s,tm)    = -inf
            c_path_sm(s,tm)       = -inf
            SMP_path_sm(s,tm)     = -inf
            risky_path_sm(s,tm)   = -inf
            safe_path_sm(s,tm)    = -inf
            share_path_sm(s,tm)   = -inf
            hcval_path_sm(s,tm)   = -inf
            else           
            
            cash_path_c(s+SS,tm)    = -inf
            c_path_c(s+SS,tm)       = -inf
            SMP_path_c(s+SS,tm)     = -inf
            risky_path_c(s+SS,tm)   = -inf
            safe_path_c(s+SS,tm)    = -inf
            share_path_c(s+SS,tm)   = -inf
            
            end if
          
             end do 
            end do
          
     
          do l=1,JJ
            cash_coh_sf(l)    = mean_noinf(cash_path_sf(:,l),SS,inf)
            cash_coh_sm(l)    = mean_noinf(cash_path_sm(:,l),SS,inf)
            cash_coh_c(l)     = mean_noinf(cash_path_c(:,l),2*SS,inf)
            
            SMP_coh_sf(l)     = mean_noinf(SMP_path_sf(:,l),SS,inf)
            SMP_coh_sm(l)     = mean_noinf(SMP_path_sm(:,l),SS,inf)
            SMP_coh_c(l)      = mean_noinf(SMP_path_c(:,l),2*SS,inf)
            
            ucshare_coh_sf(l) = mean_noinf(share_path_sf(:,l),SS,inf)
            ucshare_coh_sm(l) = mean_noinf(share_path_sm(:,l),SS,inf)
            ucshare_coh_c(l)  = mean_noinf(share_path_c(:,l),2*SS,inf)
            
            risky_coh_sf(l)   = mean_noinf(risky_path_sf(:,l),SS,inf)
            risky_coh_sm(l)   = mean_noinf(risky_path_sm(:,l),SS,inf)
            risky_coh_c(l)    = mean_noinf(risky_path_c(:,l),2*SS,inf)
         
            safe_coh_sf(l)    = mean_noinf(safe_path_sf(:,l),SS,inf)
            safe_coh_sm(l)    = mean_noinf(safe_path_sm(:,l),SS,inf)
            safe_coh_c(l)     = mean_noinf(safe_path_c(:,l),2*SS,inf)
            
            c_coh_sf(l)       = mean_noinf(c_path_sf(:,l),SS,inf)
            c_coh_sm(l)       = mean_noinf(c_path_sm(:,l),SS,inf)
            c_coh_c(l)        = mean_noinf(c_path_c(:,l),2*SS,inf)
            
            hcval_coh_sf(l)   = mean_noinf(hcval_path_sf(:,l),SS,inf)
            hcval_coh_sm(l)   = mean_noinf(hcval_path_sm(:,l),SS,inf)
        
        
			! conditional risky share 
			! single women 
            if(SMP_coh_sf(l)>0d0) then
            var=0d0
            help=0d0
            do i=1,SS
            if(SMP_path_sf(i,l)==1) var=var+1d0
            if(share_path_sf(i,l) /= -inf) help=help+ share_path_sf(i,l)
            end do
            share_coh_sf(l)   = help/var
            else
            share_coh_sf(l)   = mean_noinf(share_path_sf(:,l),SS,inf)
            end if
            
            ! single men
            if(SMP_coh_sm(l)>0d0) then
            var=0d0
            help=0d0
            do i=1,SS
            if(SMP_path_sm(i,l)==1) var=var+1d0
            if(share_path_sm(i,l) /= -inf) help=help+ share_path_sm(i,l)
            end do
            share_coh_sm(l)   = help/var
            else
            share_coh_sm(l)   = mean_noinf(share_path_sm(:,l),SS,inf)
            end if
            
            ! couples
            if(SMP_coh_c(l)>0d0) then
            var=0d0
            help=0d0
            do i=1,2*SS
            if(SMP_path_c(i,l)==1) var=var+1d0
            if(share_path_c(i,l) /= -inf) help=help+ share_path_c(i,l)
            end do
            share_coh_c(l)   = help/var
            else
            share_coh_c(l)   = mean_noinf(share_path_c(:,l),SS,inf)
            end if
        
            
          end do
          
         
         write(*,*) 'cohort variables done'
  
 
end subroutine


end module
